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Appendix A 


Preface 


The best way to acquire working knowledge on software is to practice with good 
examples. While there are numerous books and articles available, not many books are 
written for physicists who want to apply EXCEL to analyze physics problems with 
future prospects in mind. This book describes how to apply Microsoft EXCEL” to 
introductory experimental physics projects. The projects of this book are selected from 
the general physics laboratory courses at La Roche University. Measurement equip- 
ment used to acquire data are commonly found in student laboratories so that a reader 
may acquire his/her own data for analysis. This book may be used as a companion book 
for introductory laboratory courses as well as possible STEM projects. There are 
10 laboratory projects selected for this book, and each project has the following sections: 

Objectives to explain what will be learned from the project; 

EXCEL note to explain what tool and function are used for analyzing the project; 

Theory and procedure to overview the project background; 

Data analysis to show the results and the analysis method using EXCEL; 

Additional topics to describe gateways to more advanced numerical analysis. 


Topics of numerical analysis includes multiple graphs on the same EXCEL sheet, 
calculation of descriptive statistical parameters, a 3-point interpolation, the Euler 
and the Runge-Kutta methods to solve equations of motion, the Fourier transform 
to calculate the normal modes of a double pendulum, matrix calculations to solve 
coupled linear equations of a DC circuit, animation of waves and Lissajous figures, 
electric and magnetic field calculations from the Poisson equation and its 3D surface 
graphs, variational problems such as Fermat’s least traveling time principle and the 
least action principle, and drawing quantum particle trajectories using stochastic 
quantum dynamics. 

This book only uses EXCEL for all its analyses and charts, and utilizes EXCEL’s 
tools and functions useful to scientific calculations, including autofill, the interactive 
calculation option, data analysis tools, animation command, add-ins (Analysis 
ToolPak and Analysis ToolPak-VBA, and Solver), and VBA macro recording and 
editing capabilities. Each chapter is independent from each other and have the 
corresponding formulas. 

The author does not use any rigorous terminology from VBA programming in this 
book. The VBA programs created in this book are essential to future applications. 
However, just as a 5 year old boy can communicate with his limited vocabulary, these 
programs will be able to perform many scientific computations and visualizations. The 
author hopes that the readers can extend their knowledge through the examples to 
analyze everyday problems as well as more complicated physics problems. 


Summer 2019 


S I Cho 
chosl@laroche.edu 
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Book description 


A book of acquiring essential Microsoft EXCEL””s computational skills while 
analyzing introductory physics projects. Topics of numerical analysis include 
multiple graphs on the same sheet, calculation of descriptive statistical parameters, 
a 3-point interpolation, the Euler and the Runge—Kutter methods to solve equations 
of motion, the Fourier transform to calculate the normal modes of a double 
pendulum, matrix calculations to solve coupled linear equations of a DC circuit, 
animation of waves and Lissajous figures, electric and magnetic field calculations 
from the Poisson equation and its 3D surface graphs, variational calculus such as 
Fermat’s least traveling time principle, and the least action principle. Nelson’s 
stochastic quantum dynamics is also introduced to draw quantum particle 
trajectories. 

This book may be used as a companion book for introductory laboratory courses 
as well as possible STEM projects. 
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Chapter 1 


Response time of the nervous system 


1.1 Objectives 


To measure a random variable and obtain its graphical presentations including a 
histogram and a Gaussian curve: calculate the statistical variables to characterize the 
data distribution and demonstrate the central limit theorem using uniform random 
numbers. 


EXCEL note 


In chapter 1, EXCEL’s Analysis ToolPak add-in option is used to create the 
histogram. The Math&Trig functions are used to calculate the numerical values of 
statistical variables. Menus, tools, and functions used in this chapter are: 
sə Dropdown menus: [Data] > [Data Analysis] $ [Histogram]; [Chart 
Tools] o [Design]; [Formulas] > [More Functions] > [AVERAGE], 
[MODE . MULT], [MEDIAN], and [STDDEV]; and [Insert] o [Scatter 
graph]. 
% Popup menus: [Chart Tools] o [Select Data..l, and [Change 
Series Chart Type..]. 
sə Tools and Functions: Autofill, RAND (), and COUNTIF (Where do you 
want to look?, What do you want to look for?). 


1.2 Theory and procedure 


Figure 1.1 shows an experiment setup. A student determines the reaction time of his/ 
her own nervous system by catching a free-falling 30 cm ruler that is randomly 
released by another student. The catcher will catch the ruler between their index 
finger and their thumb to measure the dropped distance, x, of the ruler. This cycle is 
typically repeated 40 times per student and all of the data from several classes are 
collected for the statistical study. 


doi:10.1088/2053-2571/ab318fch1 1-1 © Morgan & Claypool Publishers 2019 


Numerical Calculation for Physics Laboratory Projects Using Microsoft EXCEL® 


Figure 1.1. Ruler dropping experiment. 


The dropped distance is given by x(t) = (1/2)gt?, where g = 9.80 ms~ is the 
gravitational acceleration. The response time of the catcher is thus calculated from 
the drop distance: t = ,/2x/g. There is some amount of uncertainty associated with 
the fluctuation of the nervous response time, which will be approximated by a 
Gaussian distribution because of the central limit theorem. 


1.3 Data analysis 
1.3.1 Histogram 


A histogram is a frequently used graphical presentation and appropriate to this 
experiment. The histogram routine appropriate to statistics is in the [Data 
Analysis] menu that is enabled with the [Analysis ToolPak] add-in option. 
Refer to appendix A1.2 for the add-in option. 

Suppose you have 2500 samples of data for the dropped distance. To make its 
histogram, open EXCEL and proceed the following steps: 

(1) In Cell A1, enter the label, Dist .x (m). 

(2) In successive cells starting with Cell A2 and moving downward, enter your 
data. This is the most time-consuming step. 

(3) In Cell B1, enter the label, Bin (cm). The label, “Bin, represents a possible 
range of data values. In this case, 1 cm to 30 cm. The numerical values of 
the bins should be 0.01 to 0.30 if the SI unit is used. 

(4) In the successive cells starting with Cell B2 and moving downward, enter 
the following bins: 1, 2, and continue to 30. Autofill can be used to enter 
the bin values. Refer to appendix A.1.1 for autofil1. 
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A B 
Dist. x (cm) Bin (cm) 


14 1 
12 2 Rows 7 to 25 
21 3 : 
= R“ are hidden for 
14 5 truncated 
11 25 : 
displa 
8 26 play 
12 27 purpose. 
8 28 
9 29 
12 30 


Data Analysis 


Analysis Tools 


Anova: Tvvo-Factor VVith Replication 
Anova: Tvvo-Factor VVithout Replication 
Correlation 

Covariance 

Descriptive Statistics 

Exponential Smoothing 

F-Test Tvvo-Sample for Variances 
Fourier Analysis 


. 
Moving Average 


Figure 1.3. [Data Analysis] window. 


(5) The initial part of the EXCEL spreadsheet should look like as shown in 
figure 1.2. Notice that, in figure 1.2, rows 7 to 25 are hidden for display 
purpose. 

(6) Click on the dropdown menu [Data] and then click on [Data Analysis] 
to display the [Analysis Tools] window (figure 1.3). Select [Histogram] 
from the dialog box. The [Histogram] window (figure 1.4) appears, and the 
cursor is now flashing in the [Input range] field. 

(7) To fill the [Input Range] field, select the data on the worksheet by 
dragging across the distance data from top to bottom (not the Bins). Drag 
the mouse to the end of data cell A2502 to highlight all data cells. The 
Input range will be entered automatically. If there are too many data 
points to highlight, the [Input Range] field can be specified as $A$2 : $A 
$2502. 

(8) Next, click on the [Bin Range] field box. To fill the [Bin Range] field, 
select the bin data on the worksheet by dragging across the data from top 
to bottom. This field is also specified as $B$2:$B$31. 

(9) Click to select the [Chart Output] option to draw a histogram. 
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Histogram Data cells ? x 
Input 
Input Range: SAS2:SAS2502] |f 
k Cancel 
Bin Range: SBS2:5BS31 m 
O tabelis Help 
Output options 


O Output Range: 
@ New Worksheet Ply: This is the default. 


O New Workbook 


mi Pareto (sorted histogram) 
İB Cumulative Percentage Don’t forget to check 
Kİ Chart Output [Chart Output] 


Figure 1.4. Preparing to draw the histogram. 


(10) Click [OK] to close the window. A histogram will be drawn on a new 
worksheet as the default is [New Worksheet Ply..]. Resize the histo- 
gram by dragging the lower edge downward. Notice that the tally (or the 
frequency) of your data is automatically calculated by EXCEL. 

(11) Click on the histogram to display the [Chart Tools] dropdown menu on 
screen. Select [Design] and then select [Add Chart Element] to enter 
[Axis Titles] and [Chart Title]. Change the name of the histogram 
appropriately (figure 1.5). 


1.3.2 Another histogram available in EXCEL 


EXCEL has another histogram routine in the [Chart] menu. Highlight the data 
cells (A2 :A2502 in this case) and follow [Insert] > [Charts]. Click on the right 
lower corner of [Charts] to display the [Insert Chart] menu (figure 1.6). Select 
the [A11 Charts] window. Click on [Histogram] and hit [OK] to obtain a 
histogram (figure 1.7). The author does not recommend this histogram because its 
“bin” values are automatically grouped and the graph is unsuitable to investigate 
statistics. 


1.3.3 Statistical variables 


There are three averages in statistics: mean, median, and mode [1]. Mean is the 
‘arithmetic’ average. Median is the ‘middle’ value in a given data set, i.e., it is the 
value separating the higher half from the lower half of the data set. Mode is the value 
that occurs ‘most often,’ or the most probable value. Although they may take the 
same numerical values in many cases, if the data distribution is asymmetric or 
skewed, they may take different values. 

Another statistical variable often used is standard deviation. It is used to indicate 
the degree of variation of a set of data. According to Chebyshev’s theorem [2], for 
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Figure 1.5. Completed histogram for the ruler-drop experiment. 
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Figure 1.6. [Insert Chart] menu. 
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Figure 1.7. Histogram from [Chart Menu] option. 


the bell-shaped Gaussian data distribution, about 68% of the data will lie within one 
standard deviation of the means, i.e., between <x>-—o and <x> + o, and, in general, 
at least 1-1/k° of the data must lie within k standard deviations on either side of the 
mean. In measurement, the standard deviation is used as the ‘error’ term. The 
measured dropped distance is then described as x(t) = <x(t)> + öx(f) where <x(t)> is 
the mean value and the error, ôx(t), is the standard deviation of the acquired data. 
Formulas for the mean value and the standard deviation are given by 


M+ + B+ ..... + Xy 
N 


Standard deviation based on sample: o = Yer —<x>). (1.2) 


Using the mean and the standard deviation, the response time of the catcher is 
given by 


t= Bisset 8) = <ft>+6, where<t>= — and ö, = 5 
\ g \ g \ 2g<x> 


As shown below, EXCEL can calculate the statistical variables using the built-in 
functions. From the data set of the histogram, we obtain the mean value <x(t)> = 12 cm 
and 6x(t) = 4.6 cm, and the dropped distance is given by x(t) = 0.12+0.05 m, and the 
reaction time is ¢ = 0.16 + 0.03 s. 


Mean: <x>= = Tax, and (1.1) 


1.3.4 Numerical calculation of statistical variables 


The calculation of mean, median, and standard deviation values is straightforward. 
On the spreadsheet created in section 3, enter Mean=, Median=, Mode=, and 
Std Dev= in Cells D1, D2, and D3, respectively in the data spreadsheet. The 
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Formulas Data Review View Help Q Te 


Ab RE Hi Ss 


Text Date & Lookup & Math& More Name 
~ Time ~ Reference ” Trig~ Functions ” | Manager 
n Library TTR statistical » İİ 


WR Engineering > 
IR Cube 2 


lo Information > 


Select this 


category. İF. Compatibility > 


İF, web , 


Function Arguments 


AVERAGE 
Number1  a2:a2502İ | = 414;1221;17;14;20;11;15;6;7:7;11;12;12;10;10;11;14; 


Number2 number 


Data Cells A2:A2502 


= 12.21671331 
Returns the average (arithmetic mean) of its arguments, which can be numbers or names, arrays, or references that contain numbers. 


Number1: numberl,number2.... are 1 to 255 numeric arguments for which you want the average. 


Formula result= 12.21671331 


Help on this function Cancel 


Figure 1.9. Calculation of average (mean). 


numerical values are to be displayed in Cells E1, E2, and E3. Excel calculates 
statistical variables using the [Statistical] functions. Follow [Formulas] > 
[More Functions] > [Statistical]. Refer to figure 1.8 for this procedure. 


Mean (figure 1.9) 

First, click on Cell E1 where the mean value is to be displayed. For calculation of the 
mean value, highlight the data points from which the mean value is calculated. Click 
on the [Formula] menu and then select [Average] in the [AutoSum] menu. 
Alternatively, follow [Formula] > [More Functions] > [Statistical] > 
[AVERAGE] to display the [Function Argument] window. Enter the data range in 
the [Number1] field. In the acquired 2500 data points, the data range is A2 :A2502. 
Hit [OK] to enter the calculated mean value in Cell E1. 
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Mode (figure 1.10) 

Begin by clicking on Cell E2 where the mode is to be displayed. For calculation of 
the mode, highlight the data points from which the mode value is calculated. Follow 
[Formula] > [More Functions] o [Statistical] — [MODE . SNGL] to display 
the [Function Argument] window. Enter the data range in the [Number1] field. 
In the 2500 data points, the data range is A2 :A2502. Hit [OK] to enter the calculated 
mode value in Cell E2. 


Median (figure 1.11) 

To display the median, click on Cell E3. For calculation of the median value, 
highlight the data points from which the median value is calculated. Follow 
[Formula] > [More Functions] o [Statistical] > [MEDIAN] to display 
the [Function Argument] window. Enter the data range in the [Number1] field. 
In the 2500 data points, the data range is A2 :A2502. Hit [OK] to enter the calculated 
median value in Cell E3. 


Function Arguments 
MODE.SNGL 
Number1 a2:a2502İ {14;12;21;17;14;20; 11;15;6;7;7; 11; 12;12;10;10;11;14; 
Number2 


Data Cells A2:A2502 


= 1 
Returns the most frequently occurring, or repetitive, value in an array or range of data. 


Number1: number1,number?2,... are 1 to 255 numbers, or names, arrays, or references that contain 
numbers for which you want the mode. 


Formula result= 13 


Help on this function 


E 
12.21671 
':a2502) l Function Arguments 


MEDIAN 


Number1 22:2250 m {14;12;21;17;14;20;11;15;6;7;7;11;12;12;10;10;11;14; 


Number2 


Returns the median, or the number in the middle of the set of given numbers. 


Number1: numberl,number2,... are 1 to 255 numbers or names, arrays, or references that contain 
numbers for which you want the median, 


Formula result= 12 


Figure 1.11. Calculation of median. 
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12.21671 
13 | Function Arguments 
12 


STDEV.S 
Std Dev= |!:a2502) i 


Numbert a2:32502) 
Number2 


Data Cells A2:A2502 


= 4,587310383 
Estimates standard deviation based on a sample (ignores logical values and text in the sample). 


Number1: numberl,number2.... are 1 to 255 numbers corresponding to a sample of a population and can 
be numbers or references that contain numbers. 


Formula result = 4.537310383 


Help on this function 


Figure 1.12. Calculation of standard deviation. 


Sample based standard deviation (figure 1.12) 

Lastly, click on Cell E4 where the standard deviation is to be displayed. For 
calculation of the standard deviation, highlight the data points from which the 
standard deviation is calculated. Follow [Formula] > [More Functions] > 
[Statistical] — [STDEV.S] to display the [Function Argument] window. 
Enter the data range in the [Number1] field. In the 2500 data points, the data range 
is A2:A2502. Hit [OK] to enter the calculated standard deviation value in Cell E4. 


1.4 Central limit theorem 
1.4.1 Gaussian curve fitting to raw data 


If there is only random error in the dropped distance, x, then its data distribution, 
fx), of a large number of data approaches the Gaussian distribution. This is a 
consequence of the central limit theorem [3]. 


exp[—(x — <x>)*/2067] 
v2zo 


where <x> is the mean and ø is the standard deviation of the data distribution: 


f(x) = (1.4) 


ae / "oaio b. Yad. (1.5) 


The histogram shown in figure 1.4 can be fitted to the Gaussian distribution with 
<x> = 12 cm and o = 4.6 cm. Table 1.1 shows the initial part of the normalized 
probability in Column C defined as [Frequency]/[Total number of events (2500 in 
this case)], the Gaussian distribution in Column D given by equation (1.4), and the 
mean and standard deviation are calculated from the raw data. Because the 
histogram of the probability is normalized, the calculated Gaussian distribution is 
normalized with a multiplication factor, m, i.e., m*f(x) where m = 0.048 in this case. 
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Table 1.1. Normalized data of the histogram and Gaussian distribution. 


A B C D 


Bin Frequency Nomalized Probility Gaussian 


1 12 0.0048 0.00436922 
2 32 0.0128 0.00728392 
3 37 0.0148 0.011577668 
4 61 0.0244 0.017545738 


Select Data Source 


Chart data range: | ='2500'!SA$2:SA$31,'2500'!SCS2:SC$31]| 


[E] Switch Row/Column 


4 Legend Entries (Series) Horizontal (Category) Axis Labels 
ES Add EZ Edit >< Remove E> Edit 


Raw data 1 


Click [Add] to add the 
Gaussian distribution chart. 


Hidden and Empty Cells Cancel 


Figure 1.13. [Select Data Source] window. 


Overlaying the two graphs on the same chart 

(12) By applying the procedure of section 3, make a histogram using Column A 
(Bin) and Column C (Normalized Probability). 

(13) Right-Click on the graph to popup the [Chart Tools] window and select 
[Select Data]. 

(14) Inthe [Select Data Source] window, click on the [Add] menu to open 
the [Edit Series] window (figure 1.13). 

(15) In the [Edit Series] window, click on the up arrow next to the [Series 
name :] box to input the graph name, e.g., Gaussian. Click on the up arrow 
next to the [Series values] box to specify the data range (figure 1.14). 

(16) Close the [Edit Series] window and the [Select Data Source] 
window (figure 1.15). Notice that the Gaussian distribution is also displayed 
as a histogram at this moment. 

(17) For a better view, change the Gaussian histogram to a continuous bell- 
shaped curve. Right-click on a data of the Gaussian distribution histogram 
to open the popup menu (figure 1.16). Select [Change Series Chart 
Type...]. 

(18) There is the [Choose the chart type and axis for your data 
series:] selection in the [Change Chart Type] window (figure 1.17). 
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Click this arrow 
to enter the 
chart name. 


Edit Series 


Click this arrow saris name 


to enter the 
data range. 


Note: ‘2500!’ is the sheet 


name that is automatically 
x 25001$D$2:5D$31) entered 


Edit Series 


Figure 1.14. [Edit Series] window. 
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Figure 1.15. [Change Series Chart Type...] window. 


Click on the selection box of the Gaussian distribution and select the 
[Line] graph. Hit [OK] to close the window. In order to see the graph 
better, it would be advised to change the graph color. Right-click on the 
Gaussian curve and then select a preferred color by clicking on the 
[Outline] menu. Figure 1.18 shows the completed Gaussian curve 
superposed on the histogram. 


1.4.2 Uniform distribution and the central limit theorem 


As another consequence of the central limit theorem, a sum of uniform random 
numbers can be approximated by a Gaussian distribution. For this demonstration, 


1-11 


Numerical Calculation for Physics Laboratory Projects Using Microsoft EXCEL® 


Click on the arrow and select a 
new chart type of the 
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Figure 1.16. Changing chart type of the Gaussian curve. 
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Figure 1.17. Selecting XY (scatter). 


EXCEL’s random number generation function, RAND (), which generates uniform 
random numbers between 0 and 1, can be used. Here the sum of uniform random 
numbers between —0.5 and +0.5 is created by consecutively generating 24 random 
numbers of RAND () -0.5: 
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Figure 1.18. Gaussian curve fitting to histogram. 


[RANDO — 0.5] + [RANDO — 0.5] + [RANDO — 0.5] 


+ [RAND() — 0.5] + [RANDO — 0.5) 

+ [RANDO — 0.5] + [RANDO — 0.5] + [RANDO — 0.5] 

+ [RANDO — 0.5] + [RANDO — 0.5] 

+ [RANDO — 0.5] + [RANDO — 0.5] + [RANDO — 0.5) ə 
+ [RAND() — 0.5] + [RANDO — 0.5] 

+[RAND( — 0.5] + [RANDO — 0.5] + [RANDO — 0.5] 

+ [RAND() — 0.5] + [RANDO — 0.5] 

+ [RANDO — 0.5] + [RANDO — 0.5] + [RANDO — 0.5] 

+ [RAND() — 0.5]. 


Uniform random number distribution (table 1.2) 
(19) Enter =RAND ( ) “RAND( ) +RAND () +RAND ( ) +RAND ( ) +RAND ( ) 
+RAND ( ) +RAND ( ) +RAND ( ) +RAND ( ) +RAND ( ) +RAND ( ) 
+RAND ( ) +RAND ( ) +RAND ( ) +RAND ( ) +RAND ( ) +RAND ( ) 
+RAND ( ) +RAND ( ) +RAND ( ) +RAND ( ) +RAND ( ) +RAND ( ) -12 
in Cell A2. 
(20) Click on the small black square at the lower right corner of Cell A2 and 
drag to Cell A2001 to generate 2000 random numbers (1.6) to each cell. 
(21) Autofill to generate numbers from —10 to +10 with an increment 0.1 in 
the Column B: enter 10 in Cell C12. Enter =B2-0.1 in Cell B3, and click 
and drag the small black square at the lower right corner of Cell B202 
generate —10 to +10 by the increment 0.1. 
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Table 1.2. Sum of 24 random numbers and Gaussian distribution. 


A B G D E F H | 


1 [24 randomlx-value Data # Freq Probab Gauss 


2 
3 
4 
5 


0.891966 -10 
-0.67092 -9.9 
-0.48086 -9.8 
3.408043 -9.7 


0 2.23542E-13 Mean= -0.0149 
0 3.72142E-13 Sigma= 1.396172 
0 6.16355E-13 
0  1.0156E-12 


0 
0 
0 
0 


(22) Column C is the data number for taking the tally of the generated 2000 random 
numbers. Autofil1 to generate the numbers from 1 to 201 as the data number. 

(23) Take a tally of the generated 2,000 random numbers in Column D. Write 
the following statement =COUNTIF ($A$2:$A$2001, "<="&B2) in Cell 
D2. This statement counts the frequency, 1.e., how many times a particular 
range of numbers, as specified the value of Cell B2 or less occurred in the 
2,000 random numbers stored in Cells A2 to A2001. That is, the frequency 
in Cell D2 is the number of the random numbers that is equal to or less 
than the value of Cell B2, 1.e., x < —10. 

(24) Next, write another statement =COUNTIF ($A$2:$A$2001, "<="&B3) - 
COUNTIF ($AŞ2:$A$2001, "<="&B2) in Cell D3. This statement means 
that the value in Cell D3 is the frequency of the random numbers equal to or 
less than the value of Cell B3 subtracted by the frequency of the random 
numbers equal to or less than the value of Cell B2, i.e., —10 < x <-9.9, 
Autofill to fill Cell D3 to Cell D201. 

(25) Write =E1/2000 in Cell E1 for calculating probability of occurrence. 
Autofill from Cell E1 to Cell E201. 

(26) Taking similar to steps (12)-(18) of section 3, a Gaussian distribution function 
can be fitted to the distribution the uniform random numbers generated just 
now. First, find the mean value and the standard deviation in appropriate 
cells, e.g., the mean value in Cell I2 and the standard deviation in Cell I3. 

(27) Write = (EXP (-( (B2-$I$2) 42) ) / (2*$1$342) ) / ($I$3*SOQRT(2*pi())) 
*0.11 in Cell F1 as the Gaussian distribution function given by expression (1.6). 
Autofill from Cell F1 to Cell F201. Notice that the Gaussian distribution is 
normalized with the multiplication factor 0.1. 

(28) Highlight Column B. While pressing the <Ctr1> key, highlight Columns 
E and F. Insert [XY scatter graphs with smooth lines] only on 
the same graph area. Right-click on the graph to popup the [Chart 
Tools] window. For better view, select (Change Chart Type] to 
change the chart type of the uniform distribution to [Scatter (data 
point only) ]. Refer to figure 1.19 for this step. 

(29) Figure 1.20 is the completed graph. Repeatedly press <F9> to generate 
different sets of random numbers and observe how the distribution fits to 
the Gaussian curve. 
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Figure 1.19. Selecting scatter graph (data point only) for Gaussian distribution. 
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Figure 1.20. Gaussian curve fitting to the uniform random number distribution. 
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Chapter 2 


Constant acceleration motion 


2.1 Objectives 


To obtain the acceleration of a glider that slides on a tilted air track, the displace- 
ment and the instantaneous velocity of the glider are measured to calculate the 
accelerations along the tilted surface. The gravitational acceleration is also calcu- 
lated from the measured acceleration. 


EXCEL note 


In chapter 2, two or more graphs of different data sets are plotted on the same chart 
area. Menus, tools and functions used in this chapter are: 
“+ Dropdown menu: [Insert] > [Scatter graph]; and 
[Chart Tools] > [Add Chart Element] > [Trendline] > [Select 
Data]. 
“+ Popup menu: [Chart Tools] > [Select Data]. 
sə Tool: Autofill. 


2.2 Theory and procedure 


Figure 2.1 shows a glider sliding down on a tilted air-track. The glider has a constant 
acceleration determined by the tilted angle and the gravitational constant. Assuming 
the glider starts at rest, the displacement, Ax, and the instantaneous velocity, v, can 
be calculated by 


Ax = zan and v = at where a is the acceleration. (2.1) 


The acceleration can be obtained from the coefficient of the t?-term of the 
[Displacement vrs time] graph, the slope of the [Displacement vrs (time)*] graph, or 


doi:10.1088/2053-2571/ab318fch2 2-1 © Morgan öz Claypool Publishers 2019 


Numerical Calculation for Physics Laboratory Projects Using Microsoft EXCEL® 


; —.—,———— 
firi 


40 
UC ly?) 


Figure 2.1. Air-track and glider system. 


Table 2.1. Data set 1 of the air-track glider. 


4 Ax(m) #1(sec) #2(sec) #3(sec) #4(sec) #5(sec) #6(sec) #7(sec) Avet tA2 Ax (m) v=2(Ax/t) 
a5 3.46 3.34 3.39 3.33 3.4 3.34 3.35” 3.372857 11.37617 1.5 
1.4 3.19 3.24 3.23 3.21 3.2 3.27 3.24 3.225714 10.40523 1.4 
1.3 3.18 3.15 3.17 3.14 3.12 3.25 3.1 3.158571 9.976573 1.3 
1.56 1.51 2.2801 
1.13 1.152857 1.32908 
0.82 0.874286 0.764376 


Table 2.2. Data set 2 of the air-track glider. 


2.22” 2.212857 4.896737 
2.19 2.148571 4.616359 
2.04 2.055714 4.225961 


0.97 1.032857 1.066794 
0.88 0.845714 0.715233 
0.3969 


the slope of the [Velocity vrs time] graph. For constant acceleration motion, the 
instantaneous velocity can be expressed using the arithmetic average velocity 


AX Vinitial + Vfinal 


Vavereage = = . 2.2 
m.” 2 C2) 


Thus, the sliding glider starting at rest, (Vinita = 0), has the instantaneous velocity, 
Ufinal = V, after sliding the distance, Ax, is given by v = 2Uaverage = 2(AX/t). 

Tables 2.1 and 2.2 list two data sets that used two different tilted angles. In these 
tables, the data in Column I show the average of the seven measurements for each 
displacement from 1.5-0.10 m. The EXCEL’s average calculation function intro- 
duced in section 1.3.4 is also applied here. Notice Rows 8-16 are hidden for 
truncated display of these tables. Also, notice Column K is identical to Column A. 
The reason for adding Column K is that, by default, EXCEL reads a left-side column 
as the x-value and a right-side column as the y-values. There is a command for 
swapping the x- and the y-values after creating a graph, but arranging the raw data 
in this way is better to avoid confusion. 
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2.3 Data analysis 


2.3.1 Displacement vrs time graph 


(1) For selecting the x-value (time: horizontal) and the y-value (displacement: 
vertical), highlight Cells 15 to 119 by dragging the mouse, and then while 
pressing the <Ctrl> key, highlight Column K5 to K19. 

(2) Click on the dropdown menu [Insert] and then select the [Scatter 
(data point only) ] graph (figure 2.2). 

(3) Click any data point on the graph created just now. The dropdown menu 
[Chart Tools] appears. Click on the [Design] menu, and the [Add 
Chart Element] menu appears on the left-side of the screen. 

(4) Select [Add Chart Element]. Then, follow [Trendline] > [More 
Trendline Option..] o [Power] to find the time dependence of the 
displacement (figure 2.3). Do not forget to select [Display equation on 
chart] in order to obtain the equation of the fitted curve. 

(5) Enter the axis titles and the graph title using the [Add Chart Element] 
window. Figure 2.4 shows the completed graph of Ax = (1/2)at?. The 
acceleration of the first data set is thus a; = 2 x 0.137 = 0.274 ms”” from 
the equation of the graph. 

Adding the second graph on the same chart area (figure 2.5). 

(6) Right-click on the graph created just now to display the popup menu [Data 
Tools]. Open [Select Data] to add another graph of the second data set 2. 

(7) Click on [Add] to add the second graph. The x-values are the averaged time 
data in Cells 126 : 140 (figure 2.6), and the y-values are the displacement in 
Cells A26 :A40. Hit [OK] to plot data points. 

(8) Repeat steps 3 and 4 to obtain the quadrature equation of the second data 
set. Figure 2.7 shows the complete graph. 
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Figure 2.2. Selecting scatter (data point only) graph. 
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Figure 2.3. Trendline option. 
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Figure 2.4. Displacement vrs time graph. 
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Figure 2.5. Adding the second graph. 


EXCEL automatically enters this as you highlight 
Cells 126 to 140. Notice Sheet 1! is the 
sheet number of the selected data cells. 
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Figure 2.6. Inputting data range to be graphed. 
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Figure 2.7. Two graphs on a single graph sheet. 


Remark: for the t?-dependence, the second-order polynomial, which is the 
quadrature function in this case, may be used for the curve fitting with the intercept 
setting 0. Although incorrect from the motion under investigation, data set 1 is fitted 
to y = 0.1285x? + 0.0112x and data set 2 to y = 0.3193x? — 0.0272x, as shown in 
figure 2.7. For this uncertainty of the curve fittings, the accelerations are determined 
from the slopes of linear graphs of [Displacement vrs (time)”] and [Velocity vrs time] 
in the following sections. 
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Figure 2.8. Distance vrs (time)”. 


2.3.2 Displacement vrs (time)? graph 


(9) For selecting the x-and y-values, highlight Cells K5 to K19 of the spread- 
sheet by dragging the mouse, and then while pressing the <Ctr1> key, 
highlight Column J5 to J19. 

(10) Click on the dropdown menu [Insert]. Then select the [Scatter] graph 
(data point only). 

(11) Click on any data point on the graph, and the dropdown menu [Chart 
Tools] appears. Click on [Design] from the [Chart Tools] menu. The 
[Add Chart Element] menu appears on the left-side. 

(12) Follow [Add Chart Element] > [Trendline] > [More Trendline 
Option...) and then select linear fitting. Do not forget to check [Display 
equation on chart] to display the equation on the graph. Because the 
initial position is zero, we can also select [Set intercept]. Its default 
value is 0. 

(13) Enter the axis titles and the graph title from [Add Chart Element]. 
Figure 2.8 shows the completed graph. The equation indicates that the 
acceleration is a; = 0.265 ms”. 

(14) Repeat steps (6) and (7) to add the graph of the second data set. Repeat 
steps (11) to (13) to obtain the equation of the second graph. The 
acceleration is a, = 0.609 ms”? form the equation of the second graph. 


2.3.3 Velocity vrs time graph 


(15) For selecting the x-and the y-values, highlight Cells 15 to 119 by dragging 
the mouse, and then while pressing the <Ctrl> key, highlight Column 
L5 to L19. 

(16) Click on [Insert] and select [Scatter] graph (data point only). 

(17) Click on any data point on the graph. Click on [Design] of the [Chart 
Tools] menu. The [Add Chart Element] menu appears on the left-side. 
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Figure 2.9. Instantaneous velocity vrs time. 


(18) Follow [Add Chart Element] > [Trendline] > [More Trendline 
Option... Select the linear fitting and do not forget to select [Display 
equation on chart] in order to display the equation on the graph. 
Because the initial velocity should be zero, we can also select [Set 
intercept]. Its default value is 0. 

(19) Enter the axis titles and the graph title from the [Add Chart Element] 
menu. Figure 2.9 shows the completed graph. From the obtained equation, 
the acceleration is calculated: a, = 0.266 ms””. 

(20) Repeat the above steps (6) and (7) to add the graph of the second data set. 
Repeat steps (15)-(18) to obtain the equation for the second graph. The 
acceleration is a2 = 0.606 ms”” from the equation of the second graph. 


Assuming the air track is friction-free, the gravitational constant, g, can be also 
obtained although the g-value is sensitive to the tilted angle measurement. From the 
linear figures 2.8 and 2.9, the accelerations are a, = 0.266 ms”? for data set 1 and 
a = 0.608 ms”” for data set 2. The sinusoidal value of the tilted angle used in the 
measurement is the ratio of the raised height, A, to the distance of the points of the 
air track contacting with the table floor and the surface raised from the floor, L, 
sind = h/L, which is 3.85/142.5 = 2.70 x 107? for data set 1 and 8.89/ 
142.5 = 6.24 x 107? for data set 2. The gravitational constant is thus given by 
a,/sin@; = 9.85 ms”? and a>/sin@) = 9.74 ms *. The obtained g-value, from the time 
measured with a stop watch, is also fairly precise. 
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Chapter 3 


Equation of motion 


3.1 Objectives 


To measure the traveling distance of a projectile, and to compare the acquired data 
with the theoretical prediction and numerically solve the equations of motion by 
applying the Euler method and the Runge-Kutta method. Projectile motions with 
and without air resistance are analyzed by the Euler method, whereas a harmonic 
oscillator is analyzed by the Runge-Kutta method. A nonlinear oscillation and a 
planetary motion are also demonstrated using the Runge-Kutter method. 


EXCEL note 


In this book, the Euler method uses autofil1, and the Runge-Kutta method uses 
visual basic for applications (VBA) macro programming. Menus, tools and 
functions used in chapter 3 are: 
sə Dropdown menu: [Insert] > [Scatter graph]; and [View] > [Macros]. 
“ Popup menu: [Chart Tools] o [Select Data..], and [Change 
Series Chart Type..]. 
sə Function: Autofill, and VBA macro. 


3.2 Theory and procedure 
3.2.1 Projectile motion 


Figure 3.1 shows a launcher that fires a bearing ball in the horizontal direction to 
produce projectile motion. The traveling distance from the firing position to the 
landing spot on the floor is recorded on a large sheet of carbon paper placed on the 
floor. The initial velocity of the ball, which is the horizontal velocity of the ball upon 
leaving the launcher, can be calculated from the gravitational potential energy of the 
initial position of the ball. It should be noticed that because the ball rolls down on 
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Figure 3.1. Launcher. 


the launcher rails, it has both rotational kinetic energy and translational kinetic 
energy when it leaves the launcher. Appendix A.2 calculates the initial velocity of the 
rolling ball. 

For the launcher and the bearing ball we used, the initial velocity, vox, is given by 
Vox = Jeh , where / is the initial height of the ball from the firing position. The 
vertical motion of the ball after being fired is free-falling and the time to hit the floor 
is given by hy = (1/2)gr”, where hp is the height of the firing position from the floor. 
The horizontal motion has a constant velocity motion with the initial velocity, vox, 
and the horizontal traveling distance from the launcher to the landing spot on the 


floor is given by X = voxt = Voy. 2ho/g = ./2hih?. 


3.3 Data analysis 


For the apparatus used, xo = 0, Ay = 0.135 m hy = 0.775 m, vp, = 1.15 m gə 


Voy = 0 m, and the calculated horizontal traveling distance is 0.457 m. The measured 
distance is 0.459 + 0.02 m, confirming the rolling motion of the ball on the launcher. 
The error term of the measured distance is estimated from the distribution of landing 
positions on the paper. 


3.4 Solving equation of motion using the Euler method 
3.4.1 The Euler method for projectile motion 


Of the ball’s motion, the horizontal position x(t + A?) at time t + At, where At can be 
calculated from the position x(t) at time ¢ and the constant velocity dx/dt = vox: 
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dx dx 
At) = — d — = vox- .1 
x(t + At) = x(t) + i t an 2 00. (3.1) 


t 


The gravitational acceleration of a projectile motion satisfies a second-order differ- 
ential equation: 


d?y 
də un (3.2) 
dt? ğ 
which can be decomposed into two first-order differential equations: 
do, dy 
—=-g and — cu. 3.3 
a ... e 


The Euler method numerically solves the differential equation (3.3) by succes- 
sively calculating the position y(t + A?) at time f + Af using the position y(t) at time t 
and the derivative dy/dt at time t: 


y(t + At) = y(t) + = - At = y(t) + v (t)At, (3.4) 


t 


and 


d 
v(t + At) € v(t) + de - At € v(t) — gAt because r =g. (3.5) 


t t 


With the initial condition of the position and the velocity, yo and vox, and the time 
interval, e.g. At = 0.01 s; the Euler method generates the position of the trajectory 
using equations (3.1), (3.4), and (3.5). Table 3.1 shows the spreadsheet of the exact 
trajectory and the numerical solution of the Euler method for a ball with the initial 
horizontal velocity 1.15 m s”. and the initial height 0.775 m. 

For constructing the calculation sheet, take the following steps: 

(1) Enter the following initial parameter values: h; in Cell C2, hz in Cell C3, vox 
in Cell C4, and the Ar in Cell C5. Cells B2 to B5 and B7 to H7 are the 
parameter names. The g-value is in Cell F2 and the theoretical horizontal 
traveling distance X is calculated in Cell F3. 

(2) Enter the initial values in Row 8: B8 to H8 (i.e., B8 :H8) 

(3) Create Column B using autofill B8:B48 as the time interval t = 0 to 0.4 s 
by 0.01 s step. 

(4) Enter =C8+E8*$C$5 in Cell C9 as xı = xo + Vo* AL, 

=D8+F8*$C$5 in Cell D9 as yı € yı + V,9* At, 
=F8+$F$2*$C$5 in Cell F9 as v, = Vv, + (dv,/dt)o* At, and 
=$C$3-4.9*B942 in Cell H9 as the exact solution y = yo — (1/2)gr”. 

(5) Cell G9 is identical to Cell C9. Highlight Cells C9 to H9 and complete the 
sheet using autofil1. The solid line of figure 3.2 shows the projectile curve 
using Euler’s method and the exact solution of equations (3.1) and (3.2). 
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Table 3.1. Calculations for projectile motion. 


A B 
| Projectile motion 


h; (m)- 0.135 g (m/s?)- -9.8 
h,(m)- 0.775 X(m)= 0.45744 
Vo, (m/s) 1.15022 
At (s)- 0.01 


y Vy Vy Exactx Exact y 
0 0.775 1.15022 0 0 0.775 
0.0115 0.775| 1.15022 -0.098 0.0115 0.77451 
0.023 0.77402 1.15022 -0.196 0.023 0.77304 
0.03451 0.77206 1.15022 -0.294 0.03451 0.77059 
0.04601 0.76912 1.15022 -0.392 0.04601 0.76716 
0.05751 0.7652 1.15022 -0.49 0.05751 0.76275 
0.06901 0.7603 1.15022 -0.588 0.06901 0.75736 
0.08052 0.75442 1.15022 -0.686 0.08052 0.75099 
0.09202 0.74756 1.15022 -0.784 0.09202 0.74364 
0.10352 0.73972 1.15022 -0.882 0.10352 0.73531 


Projectile motion 


Height (m) 
o 
A 


o1 0 0.1 0.2 0.3 0.4 0.5 


Distance (m) 
Euler method == = Exact solution 


Figure 3.2. [Height vrs Distance] graph of a ball. 


Two graphs overlaid on the same graph area 
(6) Highlight Cells C8 :D48 to take the data points by the Euler method, and 
then while pressing the <Ctr1> key, highlight column H to graph both the 
Euler solution and the exact solution. For editing the legend, right click on 
the graph to popup the [Chart Tools] window, then follow [Select 
Data Source] > [Edit Series]. Figure 3.2 is the completed graphs. 
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Table 3.2. Euler’s method applied to a falling object with air resistance. 


A B © 
Air resistance; dv/dt=b-av 
a= 


dv/dt exacty Exactv 


0 1 100 0 

0.01 0.95 100 0.00975 

99.9999 0.0195 0.9025 99.9998 0.01903 
99.9997 0.02853 0.85738 99.9996 0.02786 
99.9994 0.0371 0.81451 99.9993 0.03625 
99.999 0.04524 0.77378 99.9988 0.04424 
99.9986 0.05298 0.73509 99.9984 0.05184 
99.9981 0.06033 0.69834 99.9978 0.05906 


3.4.2 A falling object with air resistance using the Euler method 


Consider a vertically falling object with air resistance that is proportional to the 
velocity of the object. The equation of motion is given by 


m—=g-kv or Z h-a (3.6) 


where m is the mass of the object, g is the gravitational constant, k is the 
proportional constant, a = k/m, and b = g/m. Assuming the initial condition 
where the object starts at rest falling from height yo, the velocity and the position 
are given by 


p= ior ane. geen a 5 en] (3.7) 
a a a 


Let us apply the Euler method to find the numerical solution and compare it with 
the exact solution. Table 3.2 shows the EXCEL spreadsheet for the falling speed and 
the position. 

The falling speed and the position are plotted as shown in figures 3.3 and 3.4. In 
these figures, the exact solution is expressed with a broken line that is changed from 
the default solid line by taking the following steps: 

(7) In the falling speed graph (figure 3.5), right-click on a data point of the exact 
solution to display the popup window [Chart Tools] (figure 3.6). Select 
[Format Data Series..] to open the [Format Data Series] window. 

(8) Click the [Fill & Line] item and then select a broken line from [Dash 


Type]. 
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Falling speed with air resistance 


0.25 
0.2 
- 0.15 
o 
o 
Q 
o 01 
0.05 Terminal velocity 
0 
0 0.5 iL 15 2 2.5 
Time 
Euler = = Exact 
Figure 3.3. Altered graph of terminal velocity. 
Falling object with air resistance 
100.05 
100 The constant velocity 
99.95 motion after reaching 
99.9 the terminal velocity. 
m 
D 99.85 
r 99.8 
99.75 
99.7 
99.65 
99.6 


0 0.5 1 1.5 2 2.5 


Time 
Euler = = Exact 


Figure 3.4. Falling object with air resistance. 


M N | © P 


1 OF 1 


Falling Speed with Air Resistance 


~ [Æ 
— = Series “Exact” 
Fill Outline 


Delete 


. Reset to Match Style 
Right-click ona 


data point. 
Select Data... 


Select this to Add Data Labels 
change the solid Add Trendline... 
line to broken line. 1 i} Format Data Series... 


Euler — — Exact | 
o 


Figure 3.5. Falling speed graph. 
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i Format Plot Area 


Plot Area Options v 


Select [Fill & Line] No line 


Solid line 


Gradient line 
© Automatic 
Color 
Transparency H 
Width 


Compound type 


panpe Select a 


Cap type broken 


Join type 


Figure 3.6. Changing a graph style. 


3.5 Runge-Kutta method 
3.5.1 Limitation of the Euler method 


Although the Euler method is simple and easy to apply, it may not accurately trace 
the motion when the acceleration is rapidly changing in time. For example, a 
harmonic oscillator, d’x/dt?? = —4x, or v = dx/dt and doldt = —4x, which has the 
exact solution, x = cos(27) if x = 1 at t = 0. The numerical solution given by the Euler 
method gradually deviates from the exact solution due to the accumulated computa- 
tional error. Figure 3.7 shows the harmonic oscillation with the Euler method and 
the exact solution where the time increment is At = 0.1. Even if a smaller time 
interval is taken, e.g., At = 0.01, the deviation from the exact solution inevitably 
increases. 


3.5.2 Algorithm of the Runge-Kutta method 


Another approach to solve the equation of motion, called the Runge-Kutta method, 
is a better choice [1, 2]. This book only describes the algorithm of the fourth order 
Runge-Kutta method. 
First order differential equation: = = f(t, x) 
(9) Define the increment of t: 4.) = ti + h where A is the time interval and i = 0, 
1, 2, 28.5 NV. 
(10) Calculate the following K-values using given £ = t; and x = xi. 
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Harmonic oscillator by Euler's method 


Displacement 


Time 
Euler 


Exact 


Figure 3.7. Harmonic oscillation using the Euler method and the exact solution. 


K= f(t, x) 
26:55 

i iK (3.8) 
2455 


Ka f(t +h, x + hK). 


(11) Calculate the next x-value: x;,; = x; + Ax at t = fizi, where the position 
increment is a weight average of the K-values, i.e., 


Ax = h(K, + 2Ky + 2K3 + K4)/6. (3.9) 


(12) Repeat the steps (9) to (11) for a pre-determined N-value (e.g., 100) to 
obtain a table of the t-values and the x-values. 


2 
Second order differential equation: oe = g(t, x, vV) where v = Z 


(13) Separate the differential equation into two equations: 


z = g(t, x, v) and v = = = f(t, x, v). (3.10) 


(14) Define the time increment: ti}; = ti + A, where A is the time interval and 
i — 0, 1, 2, ..., N. 
(15) For solving d = g(t, x, v), calculate the following L-values using given 


t= fi, x = x, and v = 0;: 
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Li = g(t, x, v); 
2 2 2 
.11 
ro + £, x+ He, —: ə. 
2 2 2 
LA — d + ” x+ hL3, v+ its} 
(16) Calculate the next v-value: oyuq = v; + Av at t = fizi where 
Av = h(L, + 2L, + 2L3 + L4)/6. (3.12) 


(17) Solving = = f(t, x, v) is essentially the same as the first order differential 
equation described above. In this case, there is also a velocity term. 
Calculate the following K-values at given t = t; x = x; and v = uz 


K = f(t, x, v) 

h hK AK 
555. :77 xi) 

h hK: hK, 3.13 
25555 at) 
255 


(18) Calculate the next x-value: xızq = x; + Ax at t = fizi where 
Ax = h(K, + 2Ky + 2K3 + K4)/6. (3.14) 


(19) Repeat the above steps (13)-(18) for a pre-determined N-value (e.g., 100) 
obtain a table of the t-values, the x-values, and the v-values, using the time 
increment, h, and the velocity and the position increments given by 
equations (3.12) and (3.14). 


The Runge-Kutta method is much easier to implement with EXCEL’s VBA. 


Refer to appendix A1.3 for enabling the VBA capability. Remember that the VBA 
code for the Rung-Kutta method is applicable to many other second order differ- 
ential equation problems because the mathematical function, g(t, x, v), is the only 
part that needs to be changed. 


3.5.3 Harmonic oscillator 


Figure 3.8 shows the Runge-Kutta solution for the harmonic oscillation, 
d”xldt” = —4x where g(t, x, v) = —4x and f(t, x) = v = dx/dt for equation (3.10). 
The exact solution, x = cos(27), is also shown for comparison. The time interval is 
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Harmonic oscillator by Runge-Kutta method 


Displacement 


Time 
Runge-Kutta == «æ Exact 


Figure 3.8. Harmonic oscillator using the Runge-Kutta method. 


At = 0.1 in these solutions. There is much improvement from the Euler method. 
Figure 3.9 lists the VBA program of this harmonic oscillator. The graph is drawn 
using [Insert] > [Charts] steps after completing the numerical calculation. 


3.5.4 Van der Pool equation 


The Van der Pool equation is a second order differential equation that describes 
oscillation with non-linear damping: 


d?x 


x... (3.15) 


where £ > 0 [3]. This equation can describe various oscillator circuits and their noise 
[4]. Appendix A.3 shows an RLC oscillator circuit that can be reduced to the Van 
der Pool equation. 
For applying the Runge-Kutta method to the Van der Pol equation, define the 
following functions: 
do 2 dx 
— = g(t, x, v) € e(1 — x”) — x and — = f(t, x,u) —v. (3.16) 
dt dt 
The VBA codes of the function decelerations, g(t, x, v) and f(t, x, v), are listed in 
figure 3.10. The solution of the equation behaves quite differently, depending on the 
e-value as we discuss below. 


Case 1 (small e): figure 3.11 shows the solution when e = 0.09. At first, the equation is 
closed to a simple harmonic oscillation with a small amplitude because the e-term is 
small. The oscillation is not a decreasing (damping) oscillation but exponentially 
increasing when the oscillation, x, is small. However, as the oscillation becomes 
larger, the x”-term is no longer negligible and the exponential increase is suppressed 
to reach a stationary oscillation. 


Case 2 (large e): figure 3.12 shows the oscillation when £ = 10. This is type of 
oscillation called the relaxation oscillation. In this case, the x”-term of the first 
derivative cannot be negligible. When the amplitude changes as x(t) = xocoswt, the 
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Sub RungeKutta() 
" Harmonic Oscillator: d(2)y/dx(2)--4x OR dv/dt=-4x and dx/dt=v 
u Tpitial value of t: 


t=0 
T Tnit ralilyalue of 3: 
2 = ik 
i rmitdallıyalmesonuv" 
v=0 
' Time increment: 
h - 0.1 
' Iteration n (n*h = range of x, 0 to 5 by step h=0.1): 
n = 1000 
' Writing labels and initial value in cells: 
Tabelis: 
Cells(3, 1) "Initial th 
Gelisi3 2) Yunun 
Celis (ay 3) = "initial y" 
Cells(3, 4) = "Delta t" 
" Initial values and increment: 
Cells (4s 1) = ¢ 
Cells(4, 2) = x 
Cells (4, 3) =v 
Cells (4, 4) =h 
Celts (6, 2) — TtT 
Cells(6, 3) = "x" 
Cells (6, 4) = "y" 
" Runge-Kutta parameters: 
For i = 0 To f 
cells t7, 2) =t 
Cells(i + 7, 3) =x 
Cells(i +7, 4) =v 
İsi x anten, xə Nn) 
ll = git; x, ov) 
Kos a x Ga day cı s aza) 
dee Ss eae, Ak Tol Yo 0 eal a TA 9? a 
k= se thl 2, bn 9 182 // 2 ye in 912 // 2) 
110007 Joe a x 2 12070) 
kas fE th 0 Se dale bo v EAEL) 
14 gie Ge Wi He A tae anı o) 


KIEZER q 0 FRE 
(AU dı di eee cal ae eH) 6 


x 
ıl 
daxl 
+++ 
SEE 


Next i 
End Sub 


Function g(t, x, v) 
"gedv/dt 
g=-4* x 
End Function 


Function f(t; x, v) 
"fedx/dt 
fev 
End Function 


Figure 3.9. VBA code for a harmonic oscillator. 


Function gla, tr X, V} 
" Parameter a is given in the function statement 


" Choose a = 0.09 for a stationary oscillation 
" Choose a=10 for a relaxation oscillation 
" g=dv/dt 

g=-a* (x * 2-1) *v-x 


End Function 


Funetion £(t, x, v) 
" fedx/dt 
f-vy 
End Function 


Figure 3.10. Functions in the VBA function code for Van der Pol equation. 
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Stationary oscillation 


2.5 
2 İncreasing 
x 15 harmonic 
z 1 oscillation 
0.5 
€ 
g 0 
®© 05 0 0 10 bD DO 
a İ 
2 1 Suppressed | 
a -1.5 increasing oscillation | 
-2 
-2.5 
Time 


Figure 3.11. Transition to stationary oscillation (e = 0.09). 


Relaxtion oscilation 


Displacement x 
o 


Time 


Figure 3.12. Relaxation oscillation (£ — 10). 


square of the amplitude becomes x(t)” = (1/2)xo"[1 + cos(26əf)l. Therefore, when e is 
large, because of the 2w-component from the ex?-term, higher odd-frequency terms, 
3@, 5@, 7@, ..., appear in the amplitude x(t) to form a non-sinusoidal wave called 
relaxation oscillation. 


3.6 Runge-Kutta method for simultaneous ordinary differential 


equations 


3.6.1 Algorithm 


The Runge-Kutta method can be also applied to simultaneous ordinary differential 
equations of the following forms: 


dy, 
dt 


dy, 
dt 


dx 
= g(t, x, Y, Vys v) dt =f (t, X, Y.» Vys Vy) 
and 
dy 
= g(t, X, Y, Vys Vy) dt = f(t, x, Y.» Vys v). 


(3.17) 
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Steps to solve the equations are: 
(20) Define the time increment: ti}ı = ti + 4 where h is the time interval and 
i= Ü, 1, 2, ..., N. 
: d 
(21) For solving Zx = g(t, x, y, Vx, Vy) and x - g(t, X, y, Vx, Vy), calculate 
the following L-values using given t = t, x = Xi y = Vi Vx = ox, and 
Vy = Üyi: 


Ly = g(t, X, Y. Ux, Vy); 
Ly =. g(t, X, Y, Ux, Vy); 


irə h ae hL,, hLyı əə hLxı əd hLnl 
x — öy 2 5 2 3 y 2 3:5 2 . Mp 2 > 
h AL hLyı hL hLyl 
—— 2 +) 2 2 , Vy + IN: 
h AL, hLy» AL yo AL» 
məl : 2 (3.18) 
3 ef + 7’ x+ 27 2 F 2 Vy + 2 
h AL, hLy» AL yo hLy» : 
a txs 2 (yd 2 s Ux + 2 » Vy 3 P 
h 
Lis edi + 2 x + Lə, y+ hLy3, vy + ALY, Oy + ht} 
h 
Ly4 = 8, t+ >? x+ AL ,3, y+ hLy3, Ox + AL ,3, Oy + hLys : 


(22) Calculate the next v-value: 0,441 € “xi + Av, and oyiza € Vyi + Ad, at t = fizi 
where 


Avy, = (Lə + 2L,ə + 21.3 + Lxa)/6 and 


3.19 
Av, = h(Lyi + 2Lyo + 2Ly3 + Lya)/6. ( ) 


(23) Solving mi = f(t, X, Y, Ux, Oy) and a =f, (t, X, Y, Ux, Vy) are essentially 
the same as the first order differential equation described above. 
(24) Calculate the following K-values at given t = ti, x = x;, and v = vi. 
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Ka = f(t, x, Y, Vx, Vy) 
Ky, = f(t, X, Y, Uys Vy) 


Ka= s|: + ” x+ a, y+ — v. + a, T zə 

Kaa A A 2 gl. 55: 2 . r . 3 - 

a= ili + a + a, y+ — Vx + ma, V, + e) (3.20) 
—. ” x+ Ka yy — eke ma oi - 

K,4 = s + x x + hK, y + hK)ş, vx + hK, vy + ix) 

Kya= rl + . x + hKy3, y + hKy3, o, + hKy3, Vy + hK) 


(25) Calculate the next x-value: xi}ı = xi + Ax and yizı = yi + Ay at t = fiq 
where 

Ax = h(Ky, + 2Ky2 + 2.3 + Ky4)/6 and (3.21) 

Ay= h(Kyı + 2Kyo + 2Ky3 + K,4)/6. 


(26) Repeat the above steps (20) to (25) for a pre-determined N-value (e.g., 100) 
to obtain a table of the f-values, the x-values, and the v-values. 


3.6.2 Planetary motion 


The Earth orbiting around the Sun is governed by universal gravity, and the 
equation of motion is given by 


5. Pes 
dt (x? + y2)3? d ” 

and (3.22) 
do, y dy _ 


where M is the mass of the Sun and G is the universal gravitational constant. It is 
convenient to use astronomical units, AU, for the planetary motion. One astro- 
nomical unit of length, 1 AU, is the average distance, Reartn, between the Sun and 
Earth. Thus, the velocity of Earth is 2xRçara/(1-year) = 2n AU year 1, and 
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Planetary orbits 


current orbit 
vy=8 AU/Yr 3 
Escaping 


-3 
Figure 3.13. Possible Earth orbits. 


GM =v Rann = 40° AU? year”. The escape velocity of an object at distance equal 
to Earth’s orbit from the Sun is VZGM [Rearth = 24/2z & 8.89 AU year”), 

For the initial condition of Earth, x = 1, y = 0, o, = 0, and v, = 2r, assuming 
Earth starts orbiting at the position (1, 0) on the xy-plane and moving in the positive 
y-direction. For the Runge-Kutta method, the following definitions are used: 


X 


G2 + yp S(t, X, Y, Ux, Vy) = Vy 
and (3.23) 


Y 
g(t, X, y) = mü ZA X, Y, Ox, v,) = Uy. 


g(t, x, y) = -GM 


Figure 3.13 shows the calculation results of the current orbit (with vy = 2x AU 
year), elliptic orbit (if vy = 8 AU year '), and escaping orbit (if vy =9 AU year!) 
where the time interval is 0.02 and there are 1000 points used for the numerical 
calculations. 


Remark: the three graphs of figure 3.12 can be overlaid from different sheets in the 
following way: 

(27) Calculate each case using separate sheets. 

(28) Draw one graph, when e.g., the one for v, = 8. 

(29) Follow the steps (6)-(8) in chapter 2 to display the popup window [Chart 
Tools]. Then open the [Select Data] window and add two more 
graphs from different sheets. Notice that as the X and Y ranges are selected 
from different sheets, EXCEL automatically enters the sheet names. 


Figure 3.14 lists is the VBA code for the orbit calculation. 
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Sub satellite () 
“ Planetary (satellite) motion using the Rung-Kutter motion. 
* Taitial valne of t: 


t=0 
' Initial position (x and y): 
2-0 
y= 


Initial value of vx=0 and vy value determine the orbit: 
vy=1 for the current orbit, vy=8 is a comet-like orbit, and vy=9 is escaping velocity 


= 0) 

vy = 8 
7 pages eee t: Change wy value for different orbits: 
EE DE vy =2z for the current orbit, 

n = 1000 vy= 8 for a comet-like orbit, and 


" Parameter GM =4*pai*2 is the astronomical unit. 
GM 39.4761 
" Writing labels and initial value in cells: 
" Labels: 
Cells(3, 2) - 
Cells(3, 3) = 
Cells(3, 4) 
Cells(3, 5) 
Cells(3, 6) "Initial vy" 
Gelils (əra = "derta ET 
Initial Values and increment read from cells: 
Cells (4, 2) =t 
Cells (4, 3) = x 
Cells(4, 4) = y 
Cells(4, 5) vx 
Cells(4, 6) = vy 
Cells (4, 7) =h 
'Parameter names: 
Cells(10, 2) = "t" 
Cells (10, 3) xü 
Cells(10, 4) = "y" 
Celis 10:50 uv 
Cells(10, 6) = "vy" 
" Runge-Kutta paramters: 
For i = 0 Ton 
Cells(i 
Cells(i 
Cells (i 
Cells(i 
Cells(i 
1x1 = gx(GM, 
lyl = gy(GM, 
kxl = fx(t, x, yY, Vx, vy) 


vy=9 for escaping from the sun. 


tn 


x" 


kyl fy(t, x, Y, VX, vy) 
Ix? — ux (EM, E ap lal o Sac: tar oe Weal 2, yt b kyi y 2 
ly2 = gy(GM, t +h / 2, x +h * kxl / 2, y + h * kyl / 2) 
kx2 Ee Fe ff an oe ae [iva 6015770 y thA yx 1 x 2 Ny b də 
kyzı EVIE ERES 2 XARRA, yra kyl Z 2r yAn E L 2 NV e ye 
1x3 qx (GM, ou nu yaxın məs kx2:7 2, y Eh oky2/92) 
ly3 = gy(GM, t +h / 2, x + h * kx2 / 2, y + h * ky2 / 2) 
kx3 = fx(t +h / 2, x +h * kx2 / 2, y +h * ky2 / 2, vx + h * 1x2 / 2, vy + h € ly2 / 2 
ky3 = fylt +h / 2, x +h kx27/:2, yt + ky2 / 2, ve h * 1x2 / 2, vy th * 1y? / 2 
1x4 = gx(GM, t + h, x + h * kx3, y + h * ky3) 
ly4 = gy(GM, t + h, x + h * kx3, y + h * ky3) 
kx4 = fx(t + h, x + h ” kx3, y + h * ky3, vx + h * 1x3, vy + h * ly3) 
ky4 = fy(t + h, x + h * kx3, y + h * ky3, vx + h * 1x3, vy + h * 1y3) 
= in 
qes oy dz s? (ULE Gee se ük ee les) ez ab) oi E 
vy = vy +h * (lyl + 2 * ly2 + 2 * 1y3 + ly4) / 6 
x xalının ERU T rıb ae ERANS E 
y=y+h* (kyl + 2 * ky2 + 2 * ky3 + ky4) / 6 
Next i 
End Sub 


Function gx(GM, t, x, y) 
"dvx/dtegx 
ga cs qx ese 0) 
End Function 


Function gy(GM, t, x, y) 
"dvy/dt-gy 
gy = -GM ” y / ((x” 25 y 2) * 1.5) 
End Function 


Function fx(t, x, y, vx, vy) 
"vxedx/dt 
fx vx 
End Function 


Function fy(t, x, y, vx, vy) 
"vyedy/dt 
fy = yy 
End Function 


Figure 3.14. VBA Code for planetary motion. 
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In chapter 4, there are examples of coupled equations of motion: a double 
pendulum and coupled oscillators, which can be solved with the Runge-Kutta 
method. 
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Chapter 4 


Periodic motions 


4.1 Objectives 


To compare the theoretical expression of the oscillation period of a physical 
pendulum through measurement: estimate the position that give the minimum 
period using data interpolation, and analyze motions of a double-pendulum and 
coupled oscillators with the Runge-Kutta method and the Fourier transform. 


EXCEL note 


In chapter 4, an interpolation method is introduced to find a peak position and the 
peak value of a curve. The Runge-Kutta method is applied to solve a double 
pendulum and coupled oscillators. EXCEL’s Fourier transform routine is also 
introduced briefly. Menus, tools and functions used in this chapter are: 
sə Dropdown menus: [Insert] > [Scatter Graph]; and 
[Chart tools] > [Design] > [Add Chart Element] > [Trendline]. 
% Popup menus: [Chart Tools]—> [Select Data], [Add Chart 
Element], [Change Series Chart Type...], and [Format Data 
series]. 
sə Data analysis tools: Autofill, Fourier transform (option), and VBA 
macros. 


4.2 Theory and procedure 


Figure 4.1 shows a meter ruler used as a physical pendulum. The ruler has small 
holes, used as the pivot for the oscillation, at 5 cm intervals. 
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Figure 4.1. A physical pendulum. 


Assuming that the holes are negligible, the period of a physical pendulum is given by 


2 
Tam | = [ot Me x. dm) (4.1) 
Mgx Mgx g 12x 


where 7 is the moment of inertia of the ruler when the pivot position is x meters away 
from the center of mass, and M is the mass of the ruler [1]. By applying the parallel- 
axis theorem, the moment of inertia is given by J = Jp + Mx, where Ip is the inertial 
moment when the rotational axis is about the center of mass. For the one-meter ruler 
used for this project, Io = (M/12)(L* + w°; where L = 1.00 m is the length, 
w =2.5x 107” mis the width, and M is the mass of the ruler, which is not required to 
be measured. Because L >> w, Ip ~ (M/12)L’ = M/12 kg m” assuming the holes of the 
ruler have negligible contribution to the mass. 


4.3 Data analysis 


Figure 4.2 shows the theoretical values of the period, 7, for x = 0.04—0.50 m by step 
0.01 m and measured values from 0.05—0.45 m by step 0.05 m. Let us construct a 
graph with discrete points from the measurement data in addition to a graph that 
shows a smooth continuous line from the theoretical values. 

(1) Highlight all data points and insert scatter graphs of the data points. The 
obtained graphs are shown in figure 4.2. Notice the theoretical curve is also 
represented as a series of dots, which will be changed into a continuous line 
in the following step. 
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1.562958 
1.552572 
1.544163 
1.537505 
1.532406 


More Scatter Charts... 


1.62209 @Series1 @Series2 
1.630083 o 


1.63821 


Figure 4.2. Theoretical curve and acquired data. 


Note: click on a data point of the obtained chart. Follow [Chart Tools] > 
[Design]. Click on the [Add Chart Element] menu and then enter the chart and 
the axis titles. 

(2) For the theoretical period, a graph with a smooth line is better. Right-click 
on the chart to select [Change Chart Type]. Then, select [Choose the 
chart type and axis for your data series:] in the displayed 
popup menu [Change Chart Type] (figure 4.3). Select [Scatter with 
Smooth Lines] for ‘Theory’ and hit [OK]. 

(3) For changing the legend, right-click on the chart to choose [Select Data 
Source] (figure 4.4). Click on the [Edit] icon to open the [Edit Series] 
window. Change the series names and hit [OK] (figure 4.5). Figure 4.6 is the 
finalized graph with the edited legend. 
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Figure 4.3. Changing chart types. 


Select Data Source 


Chart data range: ESheet1158532:5D$74İ 


[E] switch Row/Column 


Horizontal (Category) Axis Labels 


EM Measurement 


Click on [Edit] to open 
[Edit Series] window. 


Hidden and Empty Cells 


Figure 4.4. Editing legend. 
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Edit Series 


Series name: 


x "Theory”İ z Theory 


Series X values: 

=Sheet1!SB$32:SBS78 = 0.04, 0.05, 0.... 
Series Y values: 

z Sheet11$C$32:5C$78 = 2.923643365, 2... 
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Figure 4.5. Entering graph titles. 


Physical pendulum 


T (sec) 


0 0.1 0.2 0.3 0.4 0.5 0.6 
x (m) 


@ Data Theory 


Figure 4.6. Finalized graph. 


4.4 Further investigation—minimum period of a physical pendulum 


The period of the physical pendulum has a minimum period although it is very 
subtle. Let us investigate the location of this minimum period. The following 
function can be analyzed to obtain the length, xin, called the minimum point, for 
the minimum period: 


2 
re) = (72) “Le “5 (4.2) 


where a = 1/g and $ = Ip/Mg. The function f(x) and hence the period T(x) have 
minimum values at the minimum point, x = Xmin = yla = ADİM se Lİ 4/12. 
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Table 4.1. EXCEL data around a minimum point. 


x (m) (TI2a)” (s) 
0.35 0.059 050 483 
0.30 0.057 449 839 


0.25 0.058 311 455 


Interpolation by curve fitting 
Figure 4.6 shows the shallow theoretical period, and the minimum point is around 


0.30 m. Because there are a limited number of raw data points, the measurement did 
not hit the predicted minimum point, x;in. However, we can estimate it by applying 
EXCELs curve fitting routine [Trendline]. Table 4.1 shows the selected three 
data points near the minimum point. 


(4) For applying EXCEL’s curve fitting routines (trendline), click on any data 
point on the scatter chart and then click on [Design] of the [Chart 
Tools] menu. 

(5) Click on the dropdown menu [Add Chart Element] and then select 
[Trendline]. Select [More Trendline Options...] and check 
[Polynomials]. 

(6) EXCEL offers up to sixth-order polynomial curve fitting. Checking the box 
[Display equation on Chart] displays the polynomial equation on 
screen. Figure 4.7 shows the second order polynomial fitting. 

(7) In this case, we can find the peak position and the minimum value by 
applying the quadrature curve fitting. 

(i) Select 3 data points and make a scatter graph. 

(ii) Apply trendline (Polynomial and Order 2) to find the best fitted 
quadrature equation for the 3 points around a peak by checking 
[Display equation on Chart]. Visually check the curve fitting 
(figure 4.8). The equation displayed on the chart is quadrature: 
y= ax? + bx + c where a = 0.4925, b = —0.2881, and c = 0.0996. 

(iii) Change the equation to y = a(x - ŁY + (c - z) and calculate 


Xmin = Z = 0.292 m. Because the theoretical calculation gives 


X = Xmin = JBla = Jhh/M = V1/12 = 0.289 m for the l-meter 


ruler, the estimated x,,;, is fairly accurate with the theoretical value. 


4.5 More periodic motions 
4.5.1 Double pendulum 


Figure 4.9 shows a double pendulum where two masses, mı and m, are connected 
with two rods of lengths of negligible mass, denoted as £, and f, making vertical 
angles, denoted as 6, and 65. The motion of a double pendulum is described by a set 
of coupled ordinary differential equations and is very sensitive to the initial 
condition. 
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Figure 4.7. EXCEL curve fitting (Trendline). 


Quadruture curve fitting 


0.0595 
0.059 ” 
0.0585 m y = 0.4925x? - 0.2881x + 0.0996 
0.058 ae 1 
0.0575 ac E Le 
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0.24 0.26 0.28 0.3 0.32 0.34 0.36 


Figure 4.8. Three-point quadrature curve fitting near minimum period. 
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Figure 4.9. Double pendulum. 


This book focuses on the motion with the small angle approximation, i.e., 6) ex 1 
and 62 « 1. The equation of the pendulum motion can be considered in the 
following way: 

(a) The upper mass m; is a simple pendulum with additional tension in the 
lower string, T>, i.e., 


mıliğı = -mg sin 6) + T, sin(@. — 0) where T, x mg. (4.3) 


(b) The lower mass y?ə is a simple pendulum with the effective string length of 
£0, + 66>, 1.€., 


m6, + £62) = -mg sin 62 (4.4) 


where double dots above the angle variables mean the second time derivative. For 0, 
< 1 and 6) « 1, equations (4.3) and (4.4) become 


m6; + (m + mə)g0, — məgö, -Ü and tö, + LÖ, + gö, = 0. (4.5) 


Remark: the Lagrangian approach provides the following set of equations of 
motion [2]: 


(m + m)b,0 + məbö, + (m + my) gt,0; =0 and £0; + gö, + gö, — 0. (4.6) 


Although the set of equation (4.5) appears to be different from the set (4.6), by 
adding the two equation (4.5), it becomes identical to the first equation of the set 
(4.6). Therefore, both sets, (4.5) and (4.6), are equivalent. 
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Let 9, = 4e% and 0, = A,e” in the set of equation (4.5), then we obtain the 
coupled equation: 


mho? + (mi + mə)lA, + məə?) =0 and ¢,@,A, + (o?b + g)Ar = 0. (4.7) 


Eliminating A>, we obtain İmiğiğe” + (m, + m2)(6 + &)g@? + (mi + mə)gl4) = 0. In 
order to have non-zero 4: (and A>), the following equation for œ must be satisfied: 


mühite + (m + mə)(A + &)gw@? + (m + mə)g = 0, (4.8) 


which gives 


ol Ts [oon + m)(6 + b) + Sm + m4 + b)? — 4m(m + mbh |. (4.9) 
Mity 
The solutions of equation (4.9) define the normal mode frequencies, and the 
general solutions for 6, and 62 are given by the superposition of the two normal 
mode oscillations: 


6) = Aj, cos(@ t + ö)) + 4? Cos(@2t + 62) and 


4.10 
6, = Aşı COS(@ it + ön) + An cos(@t + öz) ( ) 


where Aş and ôi (i, j = 1, 2) are determined by the initial condition. 
The Runge-Kutta method introduced in chapter 3 can be applied to analyze the 
dynamics of the coupled pendulums. For this purpose, we modify the equation (4.5) to 


6, = -e" + Mg _ Pey) and ds (z "2" Jo — 0:). (4.11) 


t m, m, mı 


Let ğı = v, and 6) = v, in the VBA program we created in chapter 3, and the 
following functions of the Runge-Kutta method are used for this example: 


g(t, X, y)= {aj (2 e -»] f(t, X, Y.» Ux, Dy) = Uy, 
1 
and (4.12) 


g(t, x; y) = İn İl F 2e) mi y), Tü. X, Y, Ox, Vy) = Uy. 
2 


Suppose £) = 2) = 1 m and m, = m = 1 kg, then œ? = [2 + J2]g. The numerical 
values of normal frequencies are fı = 0.38 Hz and fọ = 0.92 Hz with these string 
lengths and masses. Figures 4.12—4.18 show the graphs of numerical solutions with 
these parameter values, but different initial conditions. Analytical calculations of 
these oscillation patterns are time consuming, while the numerical analysis provides 
us with visual results promptly. 
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4.5.2 Fourier transform 


Let us apply the Fourier transform to the computed displacements, 9: (t) and 62 (t), 
in order to find the frequencies, o,,/2z and @/2x, and the corresponding amplitudes, 
Ağ, of equation (4.11). Because EXCEL’s Fourier transform is described in detail by 
another book [3], this book briefly explains the frequency analysis method. The fast 
Fourier transform (FFT) routine is in the [Data Analysis] menu that is enable 
with the Analysis ToolPak] add-in option. Refer to appendix A1.2 for the add- 
in option. 

For the double pendulum with the assumed values of the string lengths and 
masses, we take 4096 data points at a time interval of 0.1 s. Thus, the frequency 
resolution is given by 1/(4096 x 0.1) = 2.44 ms for the Fourier transform calculation. 
Table 4.2 shows an example of an FFT sheet used to find the normal mode 
frequency. 


Table 4.2. Fourier transform calculation sheet. 


Initial point of 61 
Time step 


t = 1 
Inlgialt Initial thdta initial theta2tnitilv1 Initial v2 Delt 
0.2 o 


Ö: with Hanning window 


FFT on Ə: 


dı 
HAL IFFTA|  Pend2 
070.30213825936551 0.30214 


FFT on Əz 


Dətə no.2 theta2 Freq step Pend 1 
o o 

0.1 0.0094799 0.0024 1.107£-08 0.151109299808434+1.15801704747513E-004i 0.15111 
0.2 0.0342593 0.0049 1.61£-07 -3.46897596253756E-005S-2.4870251020272E-0 3.5£-05 


FFT2 İFFT21 


1 
2 

3 0.3 00646905 0.0073 6.846E-07 -1.93211246559874£-005-3.379169190640176- 1.9E-05 

12 4 0.4 0.0883354 07-0.009954 0.0098 1.662E-06 -1.50331599866617E-005-4.37822796195831E4 1.5£-05 075.97206531667839£-002 0.08972 
5 0.5 0.094153 Í] -0.060301 0.0122 2.769E-06 -1.32113119686417E-005-5.40971322197426E( 1.3E-05 -7.0964E-08 -2.98731993077642E-002-2.32236473012875E-00! 0.02987 
6 0.6 0.0762762 

7 0.7 0.0361426 3 


0.0813 0.0146 3.229£-06 -1.22736023698783E-005-6.45848724S015S6E4 1.2£-05 -3.8273E-07 3.5662741 1690836E-006-6.17177844886896E-007 3.6£-06 


2.086-06 -1.17367048074929E-005-7.51987169073819E4 _1.2E-05 —7.6624E-07 5.24275869845593E-007-9.33458402267576E-007  1.16-06) 


Frequency step 


Initial point of 62 Öz with Hanning window 


FFT on ©,(f) 

(8) The starting data point for the FFT of each data set is selected at a nearly 
zero-angular displacement of each data set. In table 4.2, data points are in 
Cells C8:C4103 (theta1) and in Cells E12 : E4107 (theta2). 

(9) The data points are then pre-processed of the DC offset where the average 
value of the 4096 data points is subtracted from each data point, and the 
Hanning window, 0.5—0.5cos(21n/4095) where n = 0, 1, 2, ..., 4095, is 
multiplied to suppress the spectrum distortion. These pre-processes elim- 
inate a large DC spike in the frequency spectrum. In table 4.2, the pre- 
process data sets are named to be Pend1 (Cells G8:G4103) and Pend2 
(Cells J12 :J4107). 

(10) Follow [Data] > [Data Analysis] > [Fourier Analysis] in order 
to perform the Fourier transform, (figure 4.10). 
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Figure 4.10. Fourier analysis tool. 
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Figure 4.11. Fourier analysis parameter entry. 
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(11) Click on [Fourier Analysis] of the [Data Analysis] window to 
display its data selection menu and specify the parameters for Fourier 
transforms (figure 4.11). For calculating FFT on the data set 6), click on 
the [Input Range] box and enter $G$11:$G$4104. Check [Output 
Range :] to obtain the FFT on the same sheet. Click on the [Output 
Range] box and enter $H$11:$H$4104 to have FFT on 6, (f) 
values in Column H. Click on [OK] to display the Fourier Transform 
data in Column H. 

(12) Calculate absolute values, |©,(/)|, in Column I using the IMABS function. 

(13) For calculating the FFT on data set 02, repeat steps (1)-(13). 


The obtained frequency spectra clearly show the normal mode ratio of the double 
pendulum. Figures 4.12 and 4.13 show the dominant oscillation of the normal mode 
0.38 Hz. Figures 4.14-4.17 show both normal modes are nearly equally excited, 
resulting in complicated oscillatory patterns. Figures 4.18 and 4.19 show the 
dominant oscillation of the normal mode 0.92 Hz. 


Double pendulum (9:76, “0.2 rad at t=0) 


—— Upper mass (m1) 


Lower mass (m2) 


Figure 4.12. Double pendulum angles with an initial condition of 6, = 6, = 0.2 radians. 


Oscillation spetcrum when 6,=6, =0.2 rad at t=0 


400 


500 


300 


0 0.2 0.4 0.6 0.8 1 1.2 
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Upper mass -————- Lovver mass 


Figure 4.13. Oscillation spectrum of the double pendulum (9, = 6, = 0.2 rad at z = 0). 
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Double penldulum (4,=0.2 and 6,=0 at 0) 


—— Uppler mass (m1) 


Lower mass (m2) 


Figure 4.14. Double pendulum angles with an initial condition of 0; = 0.2 and 6, = 0 radians. 


Oscillation spectrum when 6,=0.2 and 6,=0 at t=0 


N 
wo 
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Figure 4.15. Oscillation spectrum of the double pendulum (9) = 0.2 and 6 = 0 rad at £ = 0). 


Double pendulum (6,70 and @,=0.2 at t=0 ) 


— Upper mas (m1) Lovver mass (m2) 


Figure 4.16. Double pendulum angles with an initial condition of 0; = 0 and 6, = 0.2 radians. 
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Oscillation spectrum when 6,=0 and 6,=0.2 at t=0 
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Figure 4.17. Oscillation spectrum of the double pendulum (6; = 0 and ö, = 0.2 rad at t = 0). 


Double pendulum (0,=0.2 and 6,--0.2 at t=0 ) 


—— Upper mass (m1) Lovver mass (m2) 


Figure 4.18. Double pendulum angles with an initial condition of 0; = 0.2 and 6) = —0.2 radians. 


Oscillation spectrum when 6,=0.2 and 6,=-0.2 at t=0 
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Figure 4.19. Oscillation spectrum of the double pendulum (6, = 0.2 and 6, = —0.2 rad at t = 0). 
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4.5.3 Coupled oscillators 


Another example of well-known oscillation is a coupled oscillation where two 
masses mı and m are connected with three springs as shown in figure 4.20. 

Denoting the displacements as xı and xz of the masses from their equilibrium 
positions, their equations of motion are 


2 2 
mö = —kıxı + kb% — xi) and miz = -kh — xi) + kyo. (4.13) 


a a 


ra 70 .. 
enhon 


Figure 4.20. Coupled oscillators. 


Coupled oscillators (x,=x,=1 at t=0) 


-2 
x1 x2 
Figure 4.21. Coupled oscillators with an initial condition of xı = xə = 1. 
Normal modes (x,=x,=1 at t=0) 
2 


(x1+x2)/2 


(x1-x2)/2 


Figure 4.22. Normal modes of coupled oscillators with an initial condition of xı = xə = 1. 
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Coupled oscillators (x,=1 and x,=—1 at t=0) 
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Figure 4.23. Coupled oscillators with an initial condition of xı = 1 and xə = —1. 
Normal modes (x,=1 and x,=—1 at t=0) 
2 


(x1+x2)/2 


(x1-x2)/2 


Figure 4.24. Normal modes of coupled oscillators with an initial condition of xı = 1 and xə = -1. 


For the function statements for the coupled oscillations are 


ki HE 


g(t, X, y)  -—x —(y — x) TAs X, Y, Uys Vy) = Vx 
mı mı 
and (4.14) 
ko ks 
g(t, X, y) = —(y -x) + —y f(t, X, Y, Ux, Vy) = Vy 
Mp) My) 


by taking x = x, and y = xə. Figures 4.21, 4.23, and 4.25 show the analyzed results, 
taking the equals masses, x = mə = 1 and the equal spring constants, 
kı = kə = ki = 10. Remark: in these figures, the center of the displacement of m, 
is vertically shifted to the position 3 in these graphs. 
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Coupled oscillators (x,=1 and x,=0 at t=0) 
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Figure 4.25. Coupled oscillators with an initial condition of xı = 1 and xə = 0. 
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Figure 4.26. Normal modes of coupled oscillators with an initial condition of xı = 1 and xə = 0. 


The normal modes of this system are given by 
GQ=(m+%)/2 and gq — (x — »)/2. (4.15) 


For the coupled oscillation, the normal modes can be found without the Fourier 
transform. Figures 4.22 and 4.24 show the oscillation patterns using equation (4.15), 
where only one of the two normal modes are excited. 

Figure 4.25 shows complicated periodic patterns where xı = 1 and xə = 0 initially. 
In the oscillation system, there is energy exchange between the two masses. In this 
case, the oscillations pattern of each mass becomes the same periodic pattern with 
time delay because of the energy exchange. Figure 4.26 shows the decomposed 
oscillation patterns. The decomposition clearly indicates that the normal modes are 
equally excited with the same amplitudes. 
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Chapter 5 


Lissajous figures 


5.1 Objectives 


To create and analyze Lissajous figures using EXCEL, and animate the figures on 
screen. 


EXCEL note 


In chapter 5, a two-dimensional scatter graph is used to plot Lissajous figures. The 
figures are also animated for a more realistic experience. The animation method is 
simple and applicable to wave motion and others. Menus, tools and functions used 
in this chapter are: 

sə Dropdown menu: [Insert] > [Scatter graph]. 

% Tool and functions: Autofill. 

% VBA command: DoEvents and VBA macro. 


5.2 Theory and procedure 


Lissajous figures can be used to analyze the properties of any pair of simple 
harmonic motions that are at right angles to each other [1]. They can be displayed on 
an oscilloscope as a ‘real time’ XY plotter, and two sinusoidal signals are applied to 
both the horizontal (X) and the vertical (Y) inputs. For simplicity, we assume the 
two input signals of unit amplitudes are 


X(t) = sinam; ht +n) and Y(t) = sin(2am,f,t + ny) (5.1) 


where fo is the base frequency and my, m,, ny, and n, are the integers. For Lissajous 
figures, there are two different experiments we can conduct: 
I. Two sinusoidal signals of the same amplitude and phase but different 
frequencies give rise to a stationary Lissajous figure. This mode of oscilloscope 
operation is useful in measuring the frequency ratio of the two input signals. 
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II. Two signals of the same amplitude and frequency but a different phase angle 
will appear to ‘rotate’ as the phase difference between the two signals is 
varied. This can be used as a measure for the phase angle. 


Figure 5.1 shows a Lissajous figure on an oscilloscope in the XY mode. The 
actual frequencies input to the oscilloscope are fy = 1500 Hz and f, = 500 Hz, i.e., 
fx, = 3:1. Remark: this figure is a still photo but it is actually rotating on the 
oscilloscope screen due to the phase difference between the X and Y signals. 


5.3 Lissajous figures using EXCEL 


For EXCEL’s calculation, we set fọ = 1 and A, = A, = 1 in equation (5.1), and focus 
on the frequency ratio and the relative phase. Figure 5.2 shows Lissajous figures for 
time interval of one period (T = 1/fo = 1), time increment of 0.02, and my = 3, m, = 1, 
ny = ny = 0: x(t) = sin(6r?) and y(t) = sin(277). 

The calculated points are connected to draw the corresponding Lissajous 
figure using [Scatter and with Smooth Lines] in the following steps: 

(1) Input the appropriate names in the cells highlighted (figure 5.2). 

(2) Enter 0 in Cell A9, and =A9+0.02 in Cell A10. Autofill to Cell A50 to 
create time interval (0, 1] with the time increment 0.02. 

(3) Enter =+SIN(2*PI()*$B$3*$A9+PI()*$C$3) in Cell B9. 

(4) Enter =+SIN(2*PI()*$B$4*$A9+PI()*$C$4) in Cell C9. 

(5) Highlight Cells B9 and C9 and autofill to Cells B59 and C59. 

(6) Highlight data cells (B9:C59) to be graphed. 

(7) Click on the dropdown menu [Insert] and select the [Scatter with 
Smooth Lines] graph. 

(8) Click any data point on the graph created just now to display the dropdown 
menu [Chart Tools]. Click on the [Design] menu, and the [Add Chart 
Element] menu appears on the left-side of the screen. 

(9) Enter the chart and axis titles. 


Figure 5.3 shows the Lissajous figure where the phases are n, = 0 and n, = 0.5, 
which is close to the one observed on the oscilloscope (figure 5.1). 


Figure 5.1. Lissajous figure created on an oscilloscope. 
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Lissajous Figure 


1.500 | Chart] 


1.500 1.500 
1.500 
Figure 5.2. A Lissajous figure using EXCEL. 
B te D İz O e İl: s?” İİ. | I | J 
Lissajous Figure 
Frequency Lissajous Figure 
multiplication əə 
1.500 1.500 


1.500 


Figure 5.3. A Lissajous figure with relative phase difference. 
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1.500 


1.000 


1.500 


1.500 
X(#)=sin (2x fir) and YU) “sin (2r fot) X(t)=sin (2m ft) and Y(t) =sin (6 ft +0.57) 


1.500 


X(t) A, sin (2zr fut) and (2) = A, sin(27 f+ 0.577) X(t) =sin(4rf,t) and y(r) z sin öz fyt) 


Figure 5.4. Examples of Lissajous figures. 


Figure 5.4 are EXCEL graphs of typical Lissajous figures and the corresponding 
X and Y functions: 


5.4 Animation of graphs 
5.4.1 The idea of animating EXCEL chart 


Many readers would notice that whenever you change a cell value of a data set, its 
graph reflects the change immediately. EXCEL refreshes the graph each time data in 
a cell of the chart that is changed. If this is the case, then you might expect that by 
replacing the data set one after another and executing a VBA macro, the chart 
should be refreshed with a new data set. A [For ... Next] loop can replace data sets 
quickly, and then the chart should be animated! Well, this does not work. The reason 
for this is that a VBA macro program does not return its result to the operating 
system until it completes the whole task. That is, the graph will not be refreshed 
while executing each calculation of the [For ... Next] loop. Although the following 
is a temporary setback, EXCEL does have a VBA command, DoEvents, which 
temporarily pauses the execution of the macro to refresh the screen and execute any 
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Figure 5.5. A sine curve prior to animation. 


pending events in EXCEL [2]. Notice there should be a complete graph created prior 
to running a VBA macro, which changes the data of the chart stepwise. 


5.4.2 A propagating sine-wave 


Let us animate a sine wave, y(x, t) = sin(2rx/4 — 2nf t) where 2 is the wavelength and 
fis the frequency, to make the traveling wave on a spreadsheet. For simplicity, use 
Az 0.5 and f= 1. 

(10) Draw a chart of the sine function at t = 0, i.e., y(x, 0) = sin(2vx/2). Taking 
the x-increment as 0.01 and calculate 200 data points starting at x = 0. In the 
completed chart, autofil1 is used to calculate the data points. Enter ‘n’, 
“x”, and ‘y’ in Cells A6, B6, and C6, respectively. Enter the initial values: 0 in 
Cells A7, 0 in Cell B7, and = SIN(2*pi()*B7/0.5) in Cell C7. 

(11) Enter = 1+A8 in Cell A8, = 0.01+B7 in Cell B8, and =SIN(2*pi() 
*B8/0.5) in Cell C8. Autofill to complete calculation down to Row 
206 (200 data points). 

(12) Insert a scatter chart with a smooth line. Refer to figure 5.5 for part of the 
calculation and the complete sine wave. 

(13) Now, we will animate the chart. Click on the [View] pull down menu, go 
to [Macros], and [View macros]. Create the macro program as listed in 
figure 5.6. Notice the function, DoEvenst, which returns the 200 data 
points at each time increment defined as deltime. Executing this 
program moves the sinusoidal wave to the positive direction. Notice that 
the code listed here does not carefully consider the wave speed, v = Af. The 
chart speed may be adjusted by changing deltime=j/500. Remark: 
macro recoding, which will be introduced in chapter 10, may be used 
instead of creating the VBA code while carrying out steps (1)-(3). 
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Sub sinewave () 
Dim x (200) 
“ Amplitude 
a=1 
` Wavelength 
lamda = 0.5 
“ Frequency 
f=1 
For j = 1 To 500000 
deltime = 3/500 
For i = 0 To 199 
x(i) = 0.01 * i 
Cells(7 + i, 2) = x(i) 
Cells(7 + i; 3) = a * Sin(2 * 3.14 * x(i)/lamda — 2 * 3.14 € £ *deltime) 


Next i 
DoEvents 
Next j DoEvents willrefresh the data set at 
... each step of the [FOR ... NEXT] loop. 
Figure 5.6. VBA code for sine vvave animation. 
1.5 1.5 


1.5 


-1.5 -1.5 


Figure 5.7. Snapshots of an animated Lissajous figure (xin, = 1:2). 


5.4.3 Rotating Lissajous figures 


If two independent signal generators are used to supply the x and y signals, the 
resultant Lissajous figure rotates clockwise or counterclockwise on the oscilloscope 
display depending on whether the relative phase is increasing or decreasing. 


EXCEL’s Lissajous figures can also be rotated by changing the relative phase on 
the spread sheet as you observe on an oscilloscope, if we can display a graph from an 


intermediate data set during each execution. Figure 5.7 shows a series of snapshots 
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Figure 5.8. VBA code for the Lissajous figure animation. 


of an animated Lissajous figure. In this figure, the frequency ratio is nyin, = 1 : 2. 
The reader should watch how Lissajous figures are animated with different phases. 
Figure 5.8 lists the source code, which increases the relative phase from 0 to 2x by 
1500 steps. 


References 
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Chapter 6 
Kirchhoffs law 


6.1 Objectives 


To establish and solve coupled linear equations for the electrical current of a DC 
circuit using Kirchhoff’s theorem. The coupled equation is solved by applying the 
matrix method. 


EXCEL note 


Chapter 6 introduces the matrix operation: multiplication and inversion. EXEL 
menus used in this chapter are: 
% Dropdown menu: [Formulas] > [Math&Trig] — [MINVERSE] and [MMULT]. 


6.2 Theory and procedure 


Kirchhoff’s theorem states: 
I. The algebraic sum of the electrical current at any junction is zero: Ehr = 0. 


k 
II. The algebraic sum of the potential differences across all the elements, 
including batteries, of any closed loop is zero: X V, =0. 
k 
Using these rules, a set of linear equations of the circuit currents are acquired. 
EXCEL has matrix calculation routines that can solve these coupled equations. 


6.3 Circuit under measurement 


From Kirchhoffs law, the electrical currents, J), 25, 45, and 44, of the circuit (R) = 1.0 kO, 
Ry = 470 Q, Rə = 22 KQ, Ry = 3.6 KQ, R; = 1.5 kO, Vi = 6.3 V, and Və = 1.5 V), as 
shown in figure 6.1, satisfy the following equations: 
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c 

R ——— 
> 
I, 
v, 

L 
R, 

d 


Figure 6.1. An example of a DC circuit. 


I, — h — = 0 at the junction b, 
i+ 1, - 5 = 0 at the junction e, 
RA + Rh + Rsls = V for the loop abef, (6.1) 
Rab + R4l4 — Rh — — V; for the loop 5cde, and 
RA + RL + Ryl; + Rsls =V — V, for the loop abcdef. 


Using the given resistor values (in kQ), the above set of equations of J, (in mA) 
becomes 


L-hL- 1,=0 
h+h- L=0 
h + 0.475 + 1.55 = 6.3 (6.2) 


0.471 - 2.21 = 3.614 = 1.5 
L+ 2.25 +3.64 + 1.55 = 4.8. 


6.3.1 Matrix calculation using EXCEL 


In the matrix from, the set of equations become A-I = V, and the current values can 
be obtained by finding the inverse matrix of A, 471: I = A”..V. 


1-1 -1 0 0 h 0 
0 1 0 1 -ıl İb 0 
1 0.47 0 0 15İ-İRİ-İ63İ. (6.3) 
0 0.47 -2.2 -3.6 0 L 1.5 
1 0 22 36 15] |; 4.8 


Calculation of the inverse matrix A~ of equation (6.3) 
(1) Highlight A8 : E12 where the inverse matrix is to be displayed. 
(2) Click on the dropdown menu [Formula], and then click on the 
[Math&Trig] menu. Select [MININVERSE]. Once the [Function 
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Table 6.1. Error message while calculating the inverse matrix. 


> 


A 


Arguments] window appears, highlight Cells A2:E6 to specify the 
matrix array. 

(3) While pressing both the <Ctr1> and <Shift> keys simultaneously, click 
on [OK]. The inverse matrix A” should be then displayed in Cells A8 : E12. 
However, EXCEL indicates an error for this instance! (table 6.1) 

So why is this the case? The last equation obtained from the loop abcdef 
is actually redundant. It is equivalent to the sum of the equations from the 
loops abef and bcdf. Because the electrical current remains the same across 
a battery, 7, = /5 and J; = 14, the set of equation (6.2) is reduced to 


[, — h — İş — 0 at the junction b, 
(Ri + Rs) + Rob = V for the loop abef, and (6.4) 
(R; + Rab = Rh = V for the loop bcde. 


Using the numerical values of the resistors and the batteries, the above 
equations become the following matrix equation 4:7 = V 


İ ci =} lü 0 
2.5 0.47 0 İ-İDİ“İ 63 (6.5) 
0 — 0.47 5.8 b -1.5 
1 -1 -=l d 0 
where 4 =|2.5 0.47 O|,2=|b], and V=] 63 |. (6.6) 
0 -0.47 5.8 L —1.5 
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Calculation of I = AV 

Input the name of the coefficient matrix of the linear equations, A, in Cell 

A1, the name of the inverse matrix, A*, in Cell A6 and the matrix elements 

in Cells (A2 :C4). We will calculate the inverse matrix in Cells A7:C9. 

Enter the voltage data name, V, in Cell A11 and its numerical values in 

Cells A12 :A14. Enter the current value name, I (mA), in Cell C11 and the 

current names, I=, Te, and I3= in Cells C12, C13, and C14, respec- 

tively. The initial spread sheet will be shown in table 6.2. The red- 
rectangular boxes indicate where the inverse matrix and the calculated 
current values appear. 

(5) Highlight Cells A7 :C9 where the inverse matrix is to be displayed. 

(6) Click on the dropdown menu, [Formula], and then click on the 
[Math&Trig] menu. Select [MININVERSE]. Once the [Function 
Arguments] window appears, highlight Cells A2 : C4 to specify the matrix 
array (figure 6.2). 

(7) While pressing both the <Ctrl> and <Shift> keys simultaneously, click 
on [OK]. The inverse matrix is displayed in the cell array specified by step 
(5) (table 6.3). 


(4 


— 


Table 6.2. Preparation of calculating the correct inverse matrix. 


0.47 
-0.47 


1 
2 
3 
4 
5 
6 
Ti 
8 
9 
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Function Arguments 


MINVERSE 


Array A2:C4) £| = (1-1.-1:2.5,047,0:0.-047,5.5) 


34074235095... 
Returns the inverse matrix for the matrix stored in an array. 
Array is a numeric array with an et olumns, either a 


cell range or an array constant. 


Formula result = 0.148144123 


Help on this function 


Figure 6.2. Inputting the matrix A to be inversed. 


Table 6.3. Calculation result of the inverse matrix 47. 


>I 


A? 

0.148144 0.340742 0.025542 | 
-0.788 0.3152 -0.13586 

-0.06386 _0.025542_ 0.161404 Al 


ow OND U $ 0) (ə m 


I 
L 
L 


10 

11V I (mA) 
12 | 

13 

14 

15 | 
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an Function Arguments | Inverse | ? xo: 
2 1 -1 E] MMULT — matrix A? 
3 2.5 0.47 0 Array1 A7:C9 y 7 M r p = 0798149 722602032,0,34074235095... 
4 0 -0.47 5.8 Array2 A12:A14 “ T| = {0;6.3;-1.5) 
5 \ 
4 bi Vi It 7588718,2.189554915492.., 
6A Returns the matrix product of two arrays, an array with th o age array] and columns as 
7 | 0.148144 0.340742 0.025542 array2. j 
vector V 

8 -0.788 0.3152 -0.13586 Array2 is the first array Of trener „d must have the same number 
9 | -0.06386 0.025542 0.161404 of columns as Array2 has rows. 
10 
2 I (mA) Formula result = 2.108363676 
12 oi I= 12:A14) 

i Help on this function Cancel 
iz, z 

1 
14 -1.51 I= 

-——-— " 


Figure 6.3. Preparing matrix multiplication A””V. 


(8) The product of the inverse matrix and the voltage vector, A~'V, can be 
calculated using the function, MMULT. Highlight D12:D14 where A~'V is 
to be displayed. 

(9) Click on the dropdown menu [Formula], and then click on the 
[Math&Trig] menu. Select [MMULT]. Once the [Function Arguments] 
window appears, click on the box [Array 1] and highlight Cells A8 :D11 to 
specify the inverse matrix (4477) array 1. Click on the [Array 2] box and 
highlight Cells A12: A15 to specify the second matrix (V) in Array 2 
(figure 6.3). 

(10) While pressing both the <Ctrl> and <Shift> keys simultaneously, hit 
[OK]. The calculation result of A’V is displayed in Cells D12:D15 
(table 6.4). 


Table 6.4. Calculation result of matrix multiplication A~'V. 


0.148144 0.340742 0.025542 


-0.788 0.3152 -0.13586 
-0.06386 0.025542 0.161404 


1 2.108364 i 
1 
į 2-189555 1 
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6.4 Data analysis 


In the assembled circuit of figure 6.4, measure the potential differences across each 
resistor to establish a set of linear equations for the electrical currents. Figure 6.4 
shows an EXCEL sheet containing actual data. 

The actual resistor values and battery voltages are: V; = 6.22 V, V2 = 1.53 V, 
R) = 0.993 kO, Rə = 0.459 kO, R3 = 2.16 kO, R4 = 3.61 kO, and Rs = 1.47 kO, The 
voltage measurement across each resistor is: Vr; = 2.11 V, Vr = 1.02 V, 
Vr3 = —0.195 V, Vpa = —0.326 V, and Ve: = 3.14 V. Therefore, the current of 
each resistor is İ, = VrilR = 2.12 mA, h = Vr R2 = 2.23 mA, Iz = Val 
R3 = —0.09 mA, I, = V ral R4 = —0.09 mA, and I; = V rsl Rs = 2.16 mA. 


To 1.5V battery 


Figure 6.4. An assembled circuit. 


Table 6.5. Calculation result of the current from the measurement data. 


A B 
1 Matrix from Measurement 
2A 
3 -1 
4 
5 


Here are the 
calculated 
currents. 


0.459 
-0.459 
6 INVA 
7 | 0.147366 0.346599 0.02554 


8 | -0.78981 0.321059 -0.13688 
-0.06283 0.02554 0.162421 


12 [Const V Current 

13 0 1 2.11677 
14 6.22 I 2.206416 
15 -1.53 i -0.08965 
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Using the measurement values of the resistors and batteries, equation (6.4) 


becomes the following matrix equation 4:7 = V 


1 —1 -1 h 0 
whereA = 12.46 0.459 0 L/T-İLİ, and V=] 6.22 |. 
0 —0.459 5.77 L = 1.53 


Table 6.5 shows the matrix calculation, and the result agrees with the 
measurement. 
Interested readers should try the circuit in appendix A4. The calculation 


procedure is the same. 
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Chapter / 


Equipotential surface 


7.1 Objectives 


To measure and calculate the equipotential surface of a two-dimensional electric 


dipole 


and construct its three-dimensional contour, and to obtain the numerical 


solutions of the two-dimensional Poisson’s equation. 


EXCEL note 


In chapter 7, three-dimensional surface graphs are used to depict two-dimensional 
electric potential surfaces. The two-dimensional Poisson’s equation is also solved 
numerically using EXCEL’s iteration calculation option. Because the electric field is 
a vector field, the field direction at each calculation point is represented by an arrow. 
The menus, tools and functions used in this chapter are: 


Me 
kd 


Dropdown menu: [Insert] > [Chart] > (Change Chart Type] -> [A11 
Chart] > [Surface]. 

Popup menu: [3D Rotation]; [Chart Tools] o [Format Chart 
Areal. 

Option and function: Iteration calculation and autofill 

VBA command: Shapes .AddLine (Begin x, Begin y, End x, End y). 
Line 


ActiveSheet.Shapes.EndArrovheadStyle msoArrovhead 
Stealth 


ActiveSheet .Shapes.EndArrowheadWidth = msoArrovhead 
WidthMedium 
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7.2 Measurement procedure 


Figure 7.1 is a photo showing the equipment for an equipotential measurement using 
a sheet of conducting carbon-filled paper on which the two metal pins are placed to 
simulate an electric dipole. Students locate equipotential points and then connect 
them smoothly to draw equipotential lines. The equipotential lines are drawn on 
another sheet of white paper based on the measured equipotential points. The 
electric field patterns can also be estimated because the electric field lines are 
perpendicular to the equipotential lines. 

Alternatively, EXCEL’s [Surface] chart may be used to visualize better 
equipotential surfaces, although this is very coarse. The carbon paper has 20 x 20 
cells of a 1 cm square. For the voltage measurement, we subdivide each cell to 0.2 cm 
step, and used a total of 100 x 100 cells. The voltage step is 0.25 V from 7.0 V to 0.5 V. 
The potential values at the metal pins are 6.3 V and 0 V, respectively. Record 
potential values corresponding to the nearest cells. Figure 7.2 is the measurement 
result. In figure 7.2, for a better view, the font colors are changed alternatively after 
each measurement. Notice Row 1 and Column A of the spreadsheet are used strictly 
for numbering and not data. 


Figure 7.1. Equipotential measurement apparatus. 


Figure 7.2. Acquired data of equipotential points. 
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7.3 Data analysis 
7.3.1 3D Surface graph 


A three-dimensional graph of the electric potential and its rotated graphs provides a 
fresh view of the equipotential surface. 

(1) Highlight all data points (Cell B2 to EA101) and click on the [Insert] 
menu to display the [Charts] menu. Then click on its right bottom corner 
to display the [See All Charts] window (figure 7.3). 

(2) Once the [Change Chart Type] window opens, click on [A11 Charts] 
and select [Surface] (figure 7.4). Then select [Wire Frame 3D- 
Surface]. Hit [OK] to display a wire-frame three-dimensional surface 
graph of z = (x, y) (figure 7.5). 


Page Layout Formulas Data Review View Help Q Tell me v/hat you vvani 
E E EH 3 CA (Çi) Dsmart la ul J8- rhe is 
- pə g Go V Screenshot ~ I A xə lh- Tie I ia 
PivotTable Recommended Table | Pictures Online Shapes e Recommended | A i PivotChart 
PivotTables Pictures ~ Chats ~- k~ eae N 
Tables Illustrations Charts i m Y 
S y 


Figure 7.3. Selecting all charts. 
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Figure 7.4. Change chart type. 
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Figure 7.5. Equipotential surface of an electric dipole. 


Chart Title 


== Walls 
Outline 


Delete 
ən 
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ali Change Chart Type... 
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h select Data... 


(8 3-D Rotation.. 


Select 13-D Rotation] 
to display İFormat if Format Walls... 
Chart Area]. g 


Figure 7.6. 3-D Rotation of the equipotential surface. 


(3) Right click on the graph to display a popup menu and select [3-D 
Rotation] (figure 7.6). Enter rotational angles in the [Format Chart 


Area] menu and rotate the equipotential surface and observe it from 
various angles (figure 7.7). 


7.3.2 Calculation of the 3D equipotential surface of an electric dipole 


Let us consider the electric potential of a mathematical model of an electric dipole 
where there is a pair of charges “q (q = 1) separated by distance d: 
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Format Chart Area 


Chart Options ¥ Text Options 


oO 0 B 


> Shadow 
> Glow 


> Soft Edges 


3-D Format 


3-D Rotation 


Presets 


X Rotation 


Y Rotation 


Perspective 


Right Angle Axes 


V. Autoscale 


Depth (% of base) 


Default Rotation 


Reset 


Figure 7.7. Side-menu [Format Chart Areal]. 


1 


g(x, y) = 


Jœ- d +y J+ dI +y 


(7.1) 


In this example, we use the cells as discrete points of —10 < x < +10, —10 < y < +10, 
and d = 2.5 for calculation. Create the EXCEL sheet, as shown in figure 7.8, for the 


coordinates setting. 
x-axis: 
(i) Enter -10 in Cell C2. 
(ii) Enter =C2+0.5 in Cell D2. 
Gii) Autofill to Cell AQ2. 


y-axis: 
(i) Enter -10 in Cell B3. 
(ii) Enter =B3+0.5 in Cell B4. 
Gii) Autofill to Cell B43. 
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Column (x): -10 to +10 by 0.5 step 


Row (y): -10 to +10 by 0.5 step 


1 
2 
3 
4 
5 
6 
T 
8 
9 


= 
o 


Figure 7.8. Preparing an EXCEL sheet for a dipole field. 


Dipole potential 
(i) Click on Cell C3 and enter 
=1/SQRT( ($B3-1.25) 42+C$242)1/SQRT( ($B3+1.25) *2+C 
$242). 
(ii) Autofill to Cell AQ3 (to the right). 
(iii) Continue autofill to AQ43 (downward) to complete the sheet. 


Wire-framed 3D Surface Graph to make equipotential lines 
(4) Highlight Cell B2 to AQ43. 
(5) Follow [Insert] > [Chart] > [Wire-Frame 3D Surface] to draw a 
3D equipotential surface. 
(6) Enlarge the graph for a better view. 


Figure 7.9 depicts an example of rotated equipotential lines. 
3D Rotation 
(7) Similar to step (3), right-click on the 3D surface to display a popup window 
and select [3-D Rotation]. The [Format Chart Area] menu is 
displayed on the right side of the screen. Select the [3-D Rotation] 
item and the equipotential surface can be rotated. 


Changing the vertical scale 
(8) Click on the vertical scale of the chart to emphasize the vertical scale. Then 
right-click on the vertical scale to popup a menu window. Select [Format 
Axis] to display the side menu [Format Axis] on the right side of the 
screen (figure 7.10). 
(9) Click on the [Axis Option]. If it is not displayed, click the [histogram] 
icon (the circled icon in figure 7.10) and select [Axis option]. 
(10) Select [Primary Vertical Axis]. 
(11) Change BOUND: Min (change —4.0 to —1.0) and Max (change 4.0 to 1.0). 
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Electric dipole potential lines 
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Figure 7.9. Equipotential lines of electric dipole. 
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Figure 7.10. Side menu [Axis option]. 
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Electric dipole surface 
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Figure 7.11. Expanded view of the dipole potential lines. 
Figure 7.11 shows the expanded 3D equipotential lines. 
7.4 Further investigation 


7.4.1 Two-dimensional electric potential from Poisson’s equation 


Solution of Poisson’s equation with a given boundary condition provides for two- 
dimensional electric potential w(x, y) formed by a charge distribution p(x, y) [1]. 


FAX y) . ox, y) _ P% y) 
ax? öy? Eo 


(7.2) 


We discretize the independent variables x and y: x = iAx and y = jAy where i and j 
are integers, and the electric potential function, g(x, y), is defined as (i, j). Partial 
derivatives for a continuous and smooth function are calculated as 


plx, y) eG +1,j) - (i,j) öo(x, y) Pj) - oli- Ly) 
ox Ax ” ax Ax í 
or (7.3) 
plx, y) — edi, j+ 1- (i,j) öo(x, y) _ eli, j) - (i, j- 1) 
oy = Ax ” öy ` Ax ı 


Using the forward and backward derivatives, the second derivatives are given by 


Powy) 1 (H2 etd) - (262-1) 


Ox? Ax Ax Ax (7A) 
Pox y) | 1\(eGitD-96))_ (9G) - 91-01 | 
öy? Ay Ay Ay , 
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and the Poisson’s equation discretized for numerical calculation becomes: 


h h h 
h h h Eo 


(7.5) 


From equation (7.5), the electric potential, g(i, j) at the point (i, j), is given by 


a +1 j) + 00-1) + güya D+ 0GI-) PED ve) 
4 4eo 


gli, j) = 


With equation (7.6), the two-dimensional electric potential can be calculated from 
a given charge distribution, p(i, /), and a boundary condition Let us calculate the 
dipole field from which we obtain the 3D equipotential surface of figure 7.9. Notice 
that EXCEL has an option called, Iterative Calculation, which is ‘the 
repeated recalculation of a worksheet until a specific numeric condition is met’ [2]. 
Refer to appendix A.1.4 for enabling this option. 
Here is the iteration routine for the electric potential: 
(12) Use autofil1 to write x-coordinate values from —10 to +10 in Row 1 (the 
x-coordinates) and Column A (the y-coordinates). 
(13) Enter the boundary condition: q (410, y) =g(-10,y) =@(x,+10) = 
(x,-10)=0 for -10 x x < +10 and —10 x y x +10. Enter 0 in Cells 
B2:B42 of Column B, Cells AP2 :AP42 of Column AP, Cells B2:AP2 of 
Row 2, and Cells B42:AP42 of Row 42. 
(14) Enter the charges: +1(=/"/4e 9) in Cells R23 and AB23. 
(15) Enter =(B2+C1+D3+C4) /4 in Cell C3. (figure 7.12). 
(16) Autofill C3:C42 and then autofill C3:AQ3. Autofill along the rows 
and columns. Do not worry even if the cell values are all zero for the time 
being until reaching to Column O. Once Column P is autofilled, the cell 


=(B3+C2+D3+C4)/4 


ọ(i +1, j)+ọ(i-1,j)+ọ(i, j+) +ọ(i, j-1) 


Figure 7.12. Preparing the integrative calculation. 
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Electric dipole potential surface 
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Figure 7.13. Two-dimensional electric dipole potential surface. 


values are optimally changed by the Iteration calculation option. 


For Column Q, autofill Q3:Q21 and Q23 :Q42 to avoid Cells Q22. 
(17) Highlight A1:AP42 and insert a [Wire frame 3-D surface] graph 
(figure 7.13). 


7.4.2 Two-dimensional electric field from electric potential 


The electric field can be calculated from the electric potential: E(x, y, z) = 
-V o(x, y, Z). Because the electric-field is a vector, it needs two separate tables for 


E(x, y) and E(x, y) for a two-dimensional field. For the two-dimensional potential, 
the following approximation of position derivatives can be used. 


h h 
es ” EN 
Xx, l, + = L, I, = 1, 7 — 
0... vo yə. E j oij) , PI) - Ob I | 
y 2 h h 
__ 9Ün/ 1) - ç(i,/-1) 
2h 
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7.4.3 Spreadsheet for calculating electric dipole field 


From the spreadsheet of the electric dipole potential obtained through steps (17) to 
(22), E(x, y) and E(x, y) can be calculated. Here, rows change the x-coordinates: 
Cells B2 :AP2 correspond to the coordinate (—10, —10) to (+10, —10) by step 0.5, 
and columns change the y-coordinates: Cells A2 : A42 correspond to the coordinates 
(-10, —10) to (-10, +10). 


Calculation of E(x, y) 
(18) Use autofill and enter the x-coordinates from -10 to +10 in Cells B45: 
AP45 of Row 46 by step 0.5. Similarly, enter the y-coordinates from -10 
to +10 in Cells A45:A85 of Column A. 
(19) Enter the boundary value 0 in Cells B45:AP45, B45:B84; AP45:AP84, 
and B85:AP87. 
(20) Enter =- (C3-B3) in Cell B46 and autofill to AP46. 
(21) Highlight Cells B46 :AP46 and autofil1 to Cells B84 :AP84. 
Table 7.1 shows initial part of E,. 
Calculation of Ey(x, y): 
(22) Use autofil1 and enter the x-coordinates from —10 to +10 in Cells B89: 
AP89 of Row 89 by step 0.5. Enter the y-coordinates from —10 to +10 in 
Cells A90:A130 of Column A. 
(23) Enter the boundary value 0 in Cells B90:AP90; B90:B130; AP90: 
AP130; and B130:AP130. 
(24) Enter =- (C4-C3) in Cell B91. Autofill to AP91. 
(25) Highlight Cells B91:AP91. Autofill to Cells B129:AP129. 


7.4.4 Graphical representation of the electric field 


It is easy to create a graphical representation of the electric field vector from the 
spreadsheet data obtained in steps (23)-(30). The electric field components, E,(i, j) 
and F,(i, j), calculated at the discrete position (i, j) can be used to create a vector line 
segment whose x and y values are proportional to the directional cosines, E,/E and 
E,/E, where E is the magnitude of the magnetic field at position (i, j). 


For drawing field vector arrows, the VBA command 


Shapes .AddLine (Begin x, Begin y, End x, End y) .Line 


Table 7.1. Initial part of Ex. 


-10 “3.5 -3 -8.5 -ő “7.5 7 -6.5 -6 -55 3 -4.5 


ü ü ü ü ü ü ü ü ü ü ü ü 
o ü.ütü? ü.üüü? 0.0006 0.0006 0.0005 ü.ütüz 0.0004 0.0003 ü.üüüz 4E-üs -3E-05 
0 0.0014 0.0013 0.0013 0.0012 0.001 0.0009 0.0007 0.0005 0.0003 8E-05 -2E-üd 
ü 00021 08.002 0.0019 0.0018 0.0016 0.0014 0.0012 0.0003 0.0005 0.0001 -3E-04 
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is useful [3]. Notice how for this command the spreadsheet will use a different 
definition of (continuous) coordinates, as shown in figure 7.14. Because the 
coordinates were —10 < x < +10 and —10 x y x +10 for the calculation of the 


The center of a spread - 
sheet is about (440, 250). 


1 
2 
3 
4 
5 
6 
7 
8 
9 


= 
o 


Figure 7.14. The center of a spreadsheet. 


Sub Solidfieldline() 
' Drawing E-field in sheet 2, which is cleaned at first. 
Worksheets (2) .Activate 
ActiveSheet.Shapes.SelectAll 
Selection.Delete 
"Determining the origin at around the middle of a sheet. 
x000 = 400 
y000 — 250 
"Scale factor of field line segment. 
dd = 10 Əə 3 s o məna b 
"Minimum positions : Data points may be skipped to make a less 
x00 = -10 | dense graph, e.g., 
yoo = -10 1 2 
mci- oten. —— —......l. mm, dd. N.N For I 50 to 40 Step 2 
rer j =) us) CO). a m ın i For j = 0 to 40 Step 2 
" Read Ex and Ey from sheet 1 Se E EE EEE PEE E EE, ai 
Ex = Worksheets(1).Cells(46 + j, 2 + i) 
Ey = Worksheets (1).Cells(90 + j, 2 + i) 
Expand scale of positions 
K= (0) th Wey Gr ayy © 210) 
y = (y00 + 0.5 € 2) * 20 
Calculating the field vector direction 
If the field magnitudes need to be calculated, 
use dx=dd*Ex and dy=dd*Ey with proper scale factor dd. 
E — Sar (Ex 2 By 2) 
TE E = 0 Then 


dg 0 
dy 0 
Else 
dx = dd Ex / E 
dy = dd ” Ey / E 
End If 
" Stating point of the field line segment 
x0 =x 
yo = y 
" Ending point of the field line segment 
olm pulum dz 
Yyy R OY 


' Drawing the field line segment. 

' Point (x0, yO) is shifted to (x000+x0, y000+y0) around the new origin (x000, y000). 

With ActiveSheet.Shapes.AddLine(x000 + x0, y000 + y0, x000 + x, y000 + y).Line 
.EndArrowheadStyle = msoArrowheadStealth 
.EndArrowheadWidth = msoArrowheadWidthMedium 

End With 

Next 1 
Nexti 


End Sub 


Figure 7.15. Macro code for calculating electric field lines. 


7-12 


Numerical Calculation for Physics Laboratory Projects Using Microsoft EXCEL® 


electric potential and the electric field, the coordinates should be expanded for a 
larger graph on a spreadsheet, expanding to —200 < x < +200 and —200 < y < +200 
from —10 x x < +10 and —10 x y x +10. As a rule of thumb, the center of a 
spreadsheet is about (440, 250). This point is taken as the origin of the field 
coordinates, and the field points are shifted accordingly. 

Figure 7.15 lists the macro code for drawing electric field lines in Sheet 2, and 
figure 7.16 shows the obtained graph where the directions of the E-field vectors are 
only indicated. If preferred, data points can be skipped to make the graph less 
cluttered. Notice, because the field components are already calculated, the field 
vectors with the magnitude of the E-field components can be created by dx=dd*Ex 
and dy=dd*Ey where the scale factor, dd, should be selected properly. 
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Figure 7.16. Generated electric dipole field lines. 


Remarks: 
(i) In the above macro code, there is [With ... End With] statement. 


With ActiveSheet.Shapes.AddLine (x000 + x0, y000 + YÜ, 
x000 + x, y000 + y).Line 


.EndArrowheadStyle = msoArrowheadStealth 
. EndArrowheadWidth msoArrowheadWidthMedium 
End With 
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This is equivalent to 


ActiveSheet.Shapes.AddLine (x000 + x0, y000 + y0, x000 
+ x, y000 + y).Line 
ActiveSheet .Shapes . EndArrowheadStyle msoArrowhead 


Stealth 
ActiveSheet .Shapes.EndArrowheadWidth = msoArrowhead 
WidthMedium 


(ii) For commands specifying arrow style and shape, msoArrovheadStealth 
and msoArrowheadWidthMedium, refer to Microsoft.Office.Core [4, 5]. 
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Chapter 8 


Magnetic field profile 


8.1 Objectives 


To measure and calculate the magnetic field pattern of a solenoid coil, and calculate 
a two-dimensional magnetic field through the current wires by numerically solving 
the two-dimensional Poisson’s equation. 


EXCEL note 


In chapter 8, the numerical method introduced in chapter 7 is applied to calculate 
the two-dimensional magnetic field from the Poisson’s equation for vector potential. 
Menus, tools and functions used in this chapter are: 
sə Dropdown menu: [Insert] > [Chart] > [Change Chart Type] > [A11 
Chart] > [Surface]. 
% Popup menu: [3D Rotation]; [Chart Tools] o [Format Chart 
Areal. 
sə Option and function: Iteration calculation, and autofill. 
% VBA command: Shapes .AddLine (Begin x, Begin y, End x, End y). 


Line 

ActiveSheet .Shapes.EndArrowheadStyle = msoArrowhead 
Stealth 

ActiveSheet .Shapes.EndArrowheadWidth = msoArrovhead 
WidthMedium. 


8.2 Theory and procedure 


The magnetic field produced by a single current loop of radius R placed on the 
yz-plane, as shown in figure 8.1 is given by 
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Figure 8.1. Magnetic field due to a current loop. 


HIR? 


—. 8.1 
—— 2(xo? + RY” -. 


at the point P(xo, 0, 0) on the x axis, passing through the center of the current 
loop [1]. 

The magnetic field of a solenoid coil can be constructed by superposing the single 
loop fields. The maximum magnetic field due to the single current loop is at the 
center of the loop 


m Hol R? 
Binax = WR? (8.2) 
and the ratio of the magnetic field to the maximum field is given by 
By ( RP. 
= a (8.3) 


Bemax (x0 + R3 


Equation (8.3) is the essential magnetic field pattern, which will be used for 
calculating the solenoid field. 

Figure 8.2 is the VBA code to construct the magnetic field, B/Bmax, of a solenoid 
coil that has the following geometrical shape: bore diameter = 64.00 mm; wire 
diameter = 1.50 mm; length = 150 mm (100 turns); and 4 layers of solenoid coils of 
diameters, 68.5 mm, 71.5 mm, 74.5 mm, and 77.5 mm. 

Figure 8.3 shows the calculated magnetic field profile, B/Bmax, of the solenoid 
coil. For this size of the solenoid coil, the magnetic field pattern is noticeably non- 
uniform, even inside the solenoid. 


8.3 Measurement 


Figure 8.4 shows the apparatus used to measure the magnetic field. The sinusoidal 
current is supplied to the solenoid coil to produce the time dependent magnetic field, 
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Sub Solenoid() 
Dim Xobs (500) As Double 
Dim Bobs (500) As Double 
Dim y(5) As Double 
Dim x(100) As Double 
Dim By(100) As Double 
Wire diameter 
h=1.5 
Inner bore diameter 
RO = 64 
Maximum B-field 
Bmax = 0 
Write Parameter titles 
Cells(k + 2, 2) "x (mm) " 
Cells(k + 2, 3) "B/Bmax" 
Initialization 
For k = 0 To 4 
Bobs (k) = 
Xobs(k) = 
Next k 
Superposing B-field form current loop. 
Position range is -150 mm to +300 mm. The coil extends from 0 to +150 mm. 
For k = 0 To 450 
Horsuu.ulmmuosuun 


xıb nr 
By(i) = 0 
For j = 1 To 4 
1030330505) 5) 
Hyu “Ey dıb EYS RE (eles) Sey sy a 2 o 
Next j 
Bobs(k) = Bobs(k) + By(i) 
Nexti 
If Bobs(k) >= Bmax Then Bmax = Bobs(k): "Finding maximum B-field 
Cells(k + 3; 2) = Xobsi(k) 
Next k 


For k = 0 To 450 
Bobs (k) = Bobs (k) / Bmax 
Cells(k + 3, 3) = Bobs (k 
Next k 
End Sub 


Figure 8.2. The VBA code for a solenoid coil. 


B(@) = Bosin(@t). There is a small probe coil that detects the induced magnetic field. 
The induced emf is given by 


döş 
dt 


: : Vo 
= NABow sin(@t) = Vo sin(wt), and thus By = 
NAo 


emf = N (8.4) 
where N and A are the number of turns and the cross-sectional area of the probe coil. 
Vo is the amplitude of the emf signal. If the frequency is in the RF range (<10kHz), 
the back emf is negligible. 

Although the signal to noise ratio of the observed emf may be poor because of low 
output from the signal generator connected with the solenoid coil, the magnetic field 
profile is measurable without any additional amplifier. 

Figure 8.5 shows a measurement result superposed to the theoretical field profile. 
Outside the coil, the magnetic field quickly decreases and becomes too small to 
measure. 
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B-field profile of solenoid Coil 
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Figure 8.3. Theoretical calculation of the B-field. 


Probe coil: 
ID = 22 mm and 
5 turns. 


The sensitivity 
is 5 mV/div. 


Figure 8.4. Equipment setup for magnetic-field measurement. 


S.4 Additional study 
8.4.1 Magnetic field calculated from vector potential 


Magnetic fields can be derived from vector potential: B = V x A [2]. For static 
magnetic fields, the vector potential satisfies Poisson’s equation, -V7A = J, 


where J is current density. Using the vector potential, two-dimensional magnetic 
fields can be calculated by 
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B-field profile of solenoid coil 
1.2 


The probe coil did 
not pick up the 
spreading B-field 
outside the 
solenoid. 


B/Bmax 
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Probe coil position (mm) 
Theoretical field @ Measurement 


Figure 8.5. Measured and theoretical B-field profiles. 


A- 0A, 
B(x, y) = E) and Bax, y) = -AC” (8.5) 
öy : öx 
where the vector potential satisfies 
0?A.(x, 074-(x, 
9A/(x, y) + n = —uyl,(x, y). (8.6) 


ax? öy 
Similar to the numerical solution of the Poisson’s equation (7.6) for the electric 
field, the vector potential can be numerically calculated by 


Ai 1 m 4.4 —1 r A- i, j 1 A- i, j—1 WIAi, j 
Ai, i) = (i+ 1,j)+ A:G L z(i j+ 1) + 4-(i, j 177 5 J) (8.7) 


The magnetic field given by equation (8.5) can be calculated as the average of the 
forward and backward derivatives: 


aA, 
B(x, y) == ae 
.. "m + Ax, y) = 4/x, y) ) Ax, y) = Az = Ax, mi 
~ 2 Ax Ax 
- = A(x- A 
__ Ax + Ax, y) — A(x — Ax, y) ə Bit, /) 
2Ax 
= AI) Hal) 
2h ” 
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and 
OA, 
BAX, y) = — 

oy 

2 | 2 y + Ay) — A(x y) _ A y) — A(x = m 
2 Ay Ay (8.8) 

2Ay 
m 2h i 


8.4.2 Vector potential and magnetic field due to a pair of current wire 


The following example is the calculation of the magnetic field due to a pair of 
straight current wires having opposite current directions along the z-axis. The 
positive current, +1, is placed at x = +1.25 and the negative current —1, is placed at 
x = —1.25 in the two-dimensional space, —10 x x < + 10, —10 x y x +10. 
First of all, make sure that the [Iterative Calculation] option is enabled. 
Refer to A.1.5 for enabling this option. 
Calculation of the vector potential Ai, j) 
(1) Create a spreadsheet with a boundary condition. 
x-axis: 
1) Enter -10 into Cell B1. 
ii) Enter =C2+0.5 into Cell 2. 
ii) Autofill to Cell API. 
y-axis: 
i) Enter -10 into Cell A2. 
ii) Enter =B3+0.5 into Cell A3. 
ili) Autofill to Cell A42. 
(2) Enter 0 along the boundaries B2:AP2, B2:B42, AP2:AP42, and B42: 
AP42. 
(3) Enter current value -1 in Cell AA23 (x = —1.25) and +1 in Cell Q23 
(x = 41.25). 
(4) Enter = (A2+B1+C2+B3) /4 in Cell B2. 
(5) Complete the sheet by applying autofi11. 
Table 8.1 shows the initial part of the spread sheet and figure 8.6 shows 
the 3D surface graph of the completed vector potential. 
Calculation of the magnetic fields B,(i, j) 
(6) Create coordinate positions (i, 7) and a boundary condition 
X-axis: 
1) Enter -10 into Cell B44. 
ii) Enter =C2+0.5 into Cell C44. 
ili) Autofill to Cell AP44. 
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Table 8.1. The beginning of the spreadsheet for the vector potential. 


=(B4+C3+D4+C5)/4 


E E E G l K 
potential 
-10 -9.5 -9 -8.5 -! 75 7 6.5 -B -5.5 
ü ü ü ü ü ü ü ü ü 


il 0.0007] -0.0015 -0.0022 -0.0034 -0.004 -0.0045 -0.0049 -0.0052 
oj 0. -0.0029 -0.0043 -0.0069 -0.0081 -0.0091 -0.0099 -0.0105 


1 
0.8 
0.6 
0.4 
0.2 
0 
-0.2 
-0.4 
-0.6 
-0.8 

zT 


Figure 8.6. Vector potential of a pair of opposite-current wires. 
y-axis: 
i) Enter -10 into Cell AP45. 
ii) Enter =B3+0.5 into Cell A46. 
Hi) Autofill to Cell A85. 


(7) Enter 0 along the boundaries B45:AP45, B45:B85, AP45:AP85, and 
B485:AP85. 


(8) Enter =B3-C3 in Cell C46. 
(9) Complete the B, components (B45 :AP85) by performing autofill. 
Calculation of the magnetic fields B,(i, j) 
(10) Create coordinate positions (i, j) and a boundary condition 
x-axis: 
i) Enter -10 into Cell B89. 
ii) Enter =C2+0.5 into Cell C89. 
iii) Autofill to Cell AP89. 
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y-axis: 
i) Enter -10 into Cell A90. 

ii) Enter =B3+0.5 into Cell A91. 

ii) Autofill to Cell A130. 
(11) Enter 0 along the boundaries B45:AP45, B45:B85, AP45:AP85, and 

B485:AP85. 

(12) Enter =B2-B4 in Cell C91. 
(13) Complete the B, components (B91:AP130) by autofil1. 


8.4.3 Graphical representation of the magnetic field from a pair of current wires 


Creating a graphical representation of the spreadsheet obtained above steps (6)-(13) 
is straightforward. It is very much the same as the electric field discussed in section 


Sub Bfieldsolidline ( 


" B-field is graphed in sheet 2, which is cleaned up at first. 
Worksheets (2) .Activate 
ActiveSheet.Shapes.SelectAll 
Selection.Delete 

‘Determination of the origin on a spreadsheet: 


x000 = 400 

y000 = 200 
i 5777: .—..... Data points may be skipped to make a less 
" Staring from position (-10, -10) dense graph, e.g., 

200 = -10 For I = 0 to 40 Step 2 

yoo = -10 For j = 0 to 40 Step 2 


For i = 0 To 40 
For j = 0 To 40 
' Read the field values from sheet 1. 
Bx = Worksheets (1).Cells(45 + i, 2 + j) 
By = Worksheets (1).Cells(90 + i, 2 + j) 
" Expanding the scale by factor 20. 
s= (6300 3? (Sm. ah) 00 
yax Gadel ae Ose? Spb = 20) 
t Bofield line length and direction 
iz ren (55 2 Ee By 20 


maa hem 
dx = 0 
dy = 0 
Else 
de = da Bazan 
dy = -dd * By / 8 
End If 
' The starting point corresponding to position (i,j)of the B-field line. 
x0 = x 
y =y 
ü dıra” eyliyə Porn Weng khe B- frerd Mine: 
x — xud 
yy uy oo 


Draving the field line segment. 
Point (x0, yO) is shifted to (x000+x0, y0004y0) around the new origin (x000, y000). 
With ActiveSheet.Shapes.AddLine(x000 + x0, y000 + yO, x000 + x, y000 + y).Line 


.EndArrovheadStyle = msoArrowheadStealth 
-EndArrowheadWidth = msoArrowheadWidthMedium 
End With 
Next j 
Next i 
End Sub 


Figure 8.7. VBA macro code for calculating the magnetic field from vector potential. 
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Figure 8.8. Magnetic field pattern due to a pair of opposite-current wires. 


7.4.4. The magnetic field components, B,(i, j) and By(7, j) calculated at the discrete 
position (i, j) can be used to create a magnetic field vector line segment, whose x and 
y values are proportional to B,/B and B,/B, where B is the magnitude of the 
magnetic field at position (i, j). For drawing the arrows, a VBA command Shapes. 
AddLine (Begin x, Begin y, End x, End y) .Line can be applied [3]. 

Figure 8.7 lists the VBA code for drawing magnetic field line segments with the 
arrows in Sheet 2. For the origin and the coordinate on a spreadsheet, refer to 
section 7.4.4. Figure 8.8 is the calculation result where the orientations of the B-field 
vectors are only drawn. Notice, because the field components are already calculated, 
the vector length can be given by dx=dd*Bx and dy=-dd*By, where the scale 
factor, dd, should be selected properly. For the coordinate system on the spread- 
sheet, refer to figure 7.14 of section 7.4.4. For the [With .. End With] statement, 
and the commands for specifying arrow style and shape, refer to remarks below 
figure 7.16. 
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Chapter 9 


Law of refraction 


9.1 Objective 


To measure the index of refraction by applying Snell’s law, and apply Fermat’s 
principle of least traveling time to demonstrate Snell’s law; and analyze the least 
action principle and the lowest eigenvalue of Strum—Liouville equation. 


EXCEL note 


In chapter 9, the Solver option is applied to variational calculus problems. Menus, 
tool and function used in this chapter are: 

se Dropdown menus: [Insert] > [Scatter graph]. 

“+ Add-in option: Solver. 

sə Tools and Functions: Autofill. 


9.2 Theory and procedure 


Snell’s refraction law can be obtained from Fermat’s principle which states that the 
optical path takes the shortest traveling time [1]. Figure 9.1 shows an optical path 
across a boundary between two different media with indices of refraction at nı and 
nə, respectively. The traveling time from the point A (0, 0) to point B (/, m) via point 
M (a, y) on the boundary is given by 


2 2 b2 l= 2 
a VPHP | 4030-yy 0.1) 
0 02 


where vı = c/n, and v2 = c/n are the speed of light in the given media. 
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A(0,0) 


Figure 9.1. Optical path upon refraction. 


The least traveling time should satisfy 


= — — — - — (9.2) 
dx Uya —” vwb 4 (£ — yy’ 
From the geometry of figure 9.1, the incident and the refracted angles satisfy 
ə. 757 = sin 02. (9.3) 
afa? + y? vab? + (EC - y)? 


Therefore, the condition of the least traveling time becomes 


= sin) and 


dr _ sin@ sin 6) -o0 
dx U1 U2 l 


Because vı = c/n, and v2 = c/n, which is Snell’s law: 
nı sin 0, = Mm sin 02. (9.4) 


Figure 9.2 is a typical equipment used for measuring the incident and refracted 
angles of a plastic half cylinder. The index of refraction can be found as the slope of 
the graph of [sin 6) vrs sin 62], where 6) and 62 are the incident and refracted angles. 


9.3 Data analysis 
9.3.1 Angle measurement 


Figure 9.3 shows a measurement result of the index of refraction of the plastic 
cylinder, where the [Scatter] chart and the [Trend line] are used to find the 
linear equation. The slope value of the equation of the line determines the index of 
refraction, n = 1.49. With the obtained index of refraction of the plastic cylinder, the 
least traveling time across the boundary between air and the plastic material can be 
numerically determined. 
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Refracted angles@5 "Pe 


Figure 9.2. Measurement of Snell’s law. 


Snell's Law (sin@,=n-sin6,) 


i J? 


y = 1.4865x + 0.0114 eo”! 


0 0.2 0.4 0.6 0.8 
sin, 


Figure 9.3. Measurement data of refracted angles. 


9.3.2 Fermat’s least traveling principle using EXCEL’s Solver 


EXCELs Solver add-in performs optimized calculations. With the proper setting 
parameters for a mathematics and physics problem, the Solver optimizes the 
parameters to find numerical solutions. Let us investigate the law of refraction 
using Fermat’s principle. The obtained index of refraction of the plastic cylinder will 
be used to calculate the traveling time across the boundary between air and the 
plastic material. 

Solver setting for the refracted optical path 

First of all, refer to appendix A1.2 to add the Solver capability. Figure 9.4 shows 
a spreadsheet setting for finding the optical path from the start point (0, 0) to the end 
point (200, 200) across the boundary, at x = 100 on the boundary between air and 
another medium of n = 1.49. Notice that on the spreadsheet Columns I and J are 
legends explaining the equations for calculating the path lengths and the 
traveling time. 
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Home Insert’ Page Layout. Formulas Data | Review View 


Fe Tə A 4 à a (Bconnections 4 ( = Və ite vi “Au gim” 
in te tü ie) n ə : = Ə 9 | 


From From From r Esisting | Refresh ğ Text to Remov onsolidate What! Group Ungroup Subtotal 
Access Web Text ces- Connections All G Advanced | Columns Duplicate: n Analysis “ 7 - 


Figure 9.4. Solver option. 


As a possible pre-set or trial path, a straight line from the starting point (0, 0) to 
the ending point (200, 200) is passed through the point (100, 100) on the boundary. 
By adjusting the y-coordinate of the optical path on the boundary, the least traveling 
time will be determined. This adjustment is performed by the Solver. Notice that 
both the starting and the ending points are fixed. Once EXCEL completes 
the optimization, the incident and the refracted angles, calculated from the path 
of the starting point to the ending point through the optimized point on the 
boundary, are examined to see if they satisfy Snell’s law. 

The y-coordinate on the boundary is given in Cell D8 and the total traveling time 
is calculated in Cell G10. The speed of light in the free space is assumed to be c = 1, 
then the speed of light in air and in the medium are vı = c/n, © c and v2 = clm, 
respectively, where nı = 1.00 and m = 1.49. 

In table 9.1, the ‘trial’ optical path is a straight line drawn in red by using the 
[Scatter with Straight Lines] of the [Chart] menu to connect the 3 data 
points, Cells (B11, C11), (B12, C12), and (B13, C13). 


Table 9.1. Incident and refracted angles before optimization. 


QO “0 OY VT 5 Q) ə A 


Least time traveling 
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(1) Highlight Cells B11 to C13. 
(2) Go to the [Insert] menu and select the [Scatter] graph. 
(3) Select the [Scatter with Straight Lines] option. 


On the sheet, click on the [Data] menu and then select the [Solver] option. 
In the [Solver] window: 

(4) Enter the cell address of the total travelling time, $G$10, in the [Set 
Objective:], and the cell address of the y-coordinate on the boundary, 
$D$8, in the [By Changing Variable Cells:l. In this experiment, the 
optimization is to minimize the total traveling time. Check [Min] at the 
bottom of the [To:] option. 

(5) Default settings for other options are used. Hit [Solve] to start the 
optimization (figure 9.5). 

(6) Once the optimization window appears, hit [OK]. Notice that the default 
setting is to check [Keep Solver Solution] (figure 9.6). 

(7) The spreadsheet displays the least traveling path (table 9.2). 

(8) The angle calculation using the numerical values of Cells B16-B18 in 
table 9.2 confirms Snell’s law. 


Set Objective: 
To: 


By Changing Variable Cells: 2 . 
Total traveling time 


to be minimized 
(Cell G10). 


Subject to the © 


The y-coordinate at 
the boundary (Cell D8) 


İV) Make Unconstrained Variables Non-Negative 


Select a Solving Method: GRG Nonlinear iz) 


Solving Method 


Select the GRG Nonlinear engine for Solver Problems that are smooth nonlinear. Select the LP Simplex engine 
for linear Solver Problems, and select the Evolutionary engine for Solver problems that are non-smooth. 


Figure 9.5. Solver parameters. 
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Solver found a solution. All Constraints and optimality 
conditions are satisfied. Reports 
Answer 
(6) Keep Solver Solution This is the Sensitivity 
Limits 


O Restore Original Values default. 


C Return to Solver Parameters Dialog O Outline Reports 


Save Scenario... 


Solver found a solution. All Constraints and optimality conditions are 
satisfied. 

When the GRG engine is used, Solver has found at least a local optimal solution. 
When Simplex LP is used, this means Solver has found a global optimal solution. 


Figure 9.6. Generating Solver results. 


Table 9.2. Incident and refracted angles after optimization. 


G 


— index of referacti 


peed of light in vacuum = 1icz 


After opti 


Before optimization 


Numerical Calculation for Physics Laboratory Projects Using Microsoft EXCEL® 


9.4 Projectile motion based on the least action principle 
9.4.1 Lagrangian approach 


Other than Fermat’s principle, there are many physics laws that are expressed by the 
principle of variation where EXCEL’s Solver may be applied to find the optimized 
solution. One such example is Hamilton’s least action principle in analytical 
mechanics, which states that the condition is that a particle moves in a potential 
V(X) between two positions, xı at fı and xə at fə, in such a way that the integral, 
called the action 


t2 
s-f L(£, ddt, (9.5) 
ül 

takes the least possible value [2]. The integrant of equation (9.5) is called the 
Lagrangian function for the system. For a particle in a potential V(X), its 
Lagrangian is given by 


L(£, 0) = sini — V(%). (9.6) 


9.4.2 Projectile motion 


Let us find the trajectory of a projectile motion, assuming that a particle of mass 
10 kg is fired from the origin at the initial velocities vp, = 5.0 m s | and Voy = 8.0m s? 
in the horizontal and vertical directions, respectively. The Lagrangian is 


L(X, 0) = sino? + v,) — mgy — ym) — mgy. (9.7) 

Assume the position of the particle are measured as (10.0, —3.6) at £ — 2.0 s. The 
problem is to find the trajectory from the initial point (0, 0) to the final point 
(10.0, -—3.6) by minimizing the action (9.6). First, consider the straight line 
connecting these two positions, y(x) = 0.36x, as the “trial” function, and let 
the Solver minimize the action, starting from this trial function. Table 9.3 shows 
the EXCEL sheet that calculates the Lagrangian function and the action: 


1 (x; — Xi-1)? + y= yar 2 
L;= 5m m — mgy, and S — I Ldt x È Lidt, (9.8) 


where At = 0.1 s and the exact positions (x, y) are calculated using the exact solution: 
x(t) = 500xt = St and y(t) = Voyt — (1/2)gt* = 8t — 4.9t°. For the truncated display 
purpose, data rows 8—21 are hidden in table 9.3. 

Because an equal time increment is taken (At = 0.1 s), minimizing the action can 
be carried out by minimizing the sum of the Lagrangian functions, pes which is 


i 
calculated in Cell 125. Figure 9.7 shows the exact trajectory and the straight line as 
a trial function. 
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Table 9.3. Minimization of the Lagrangian of a projectile motion. 


10 Initial vx0= 
Velocity vyo= 


L=K-U 

9.43 444.6245 0 444.6245 
0.531413 5.314132 141.2 -17.64 158.84 
0.531413 5.314132 141.2 -35.28 176.48 
0.531413 5.314132 141.2 -52.92 194.12 
0.531413 5.314132 -317.52 458.72 
0.531413 5.314132 -335.16 476.36 
0.531413 5.314132 -352.8 494 


Total S 6973.025 


Projectile Motion 


y (m) 


12 


x (m) 


Exact Trial 


Figure 9.7. Trial function (straight line) and exact solution. 


The numerical value of XL i in Cell 125 will be minimized by changing the y’-values 


in Column D (Cells D5 :D23) while the positions (0, 0) at t = 0 s and (10.0, —3.6) at 
t = 2 s are fixed. 
The minimization steps are straightforward. In the [Solver] window: 

(1) Enter the cell address of the total travelling time, $1$25, in the [Set 
Objective:], and the cell address of the y-coordinates, $D$5:$D$23, in 
the [By Changing Variable Cells:l. The optimization is to minimize 
the sum of Lagrangian functions. Check [Min] bottom of the [To :] option. 

(2) There is no constraint, and the [Subject to the Constraints] box is 
left blank. 
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Solver Parameters 


ölər 
2 “ 


x 
Set Objective: t sis25 i 


- z 


To: O Max (9) Min O value Of: This is the sum of 


_ — By-Chənging Variable Cells: the Lagrangians to 
.— İşpss:5D$23 27 be minimized. 


Subject to the Constraints: 
Add 
Change the y’- Change 
values of the Delete 


trial function. 
Reset All 


Load/Save 


LT Make Unconstrained Variables Non-Negative 


Select a Solving GRG Nonlinear Options 
Method: 


Solving Method 
Select the GRG Nonlinear engine for Solver Problems that are smooth nonlinear. Select the LP Simplex 


engine for linear Solver Problems, and select the Evolutionary engine for Solver problems that are 
non-smooth. 


Figure 9.8. Solver parameter. 


(3) Default settings for other options are also used. Hit [Solve] to optimize the 
total traveling time. 


Figure 9.9 shows the trajectory after minimizing the Lagrangian sum and the 
exact trajectory. In the figure, the exact trajectory is expressed by a set of discrete 
points. Overall, the optimized numerical solution is accurate. Notice that through 
the choice of a better trial function it can improve the calculation efficiency for a 
better solution. 


9.5 Eigen value problems using Solver 

9.5.1 Eigenvalues of Strum—Liouville equation 

Another example of a variational problem is to find the smallest eigen values of the 
Strum-Liouville equation [3], 

du(x) 


fy] — q(x)u(x) + Ap(x)u(x) = 0 (a € x € b), (9.9) 
dx dx 
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Projectile Motion 


O e N WW 68 


12 


Vertical Position 


1 1 ü 1 
w N bəni 


Horizontal Position 


Solver 6 Theory 


Figure 9.9. Optimized path from Solver calculation and exact solution. 


which is the Euler-Lagrange equation for the functional 


b 2 
T[u(x)] = / (=) .— (9.10) 


Alternatively, vve may minimize 

Tus =f HEC (9.11) 
subject to the condition 

J[u(x)] = I . pu?dx — constant. (9.12) 


Then, the eigen value 2 enters as a Lagrange multiplier. This variational procedure is 
equivalent to minimizing 


b 
2 
K[u(x)] = 1lu(x)) - / (p(duldx)* + qu2)dx 
J[u(x)] oe 


and the stationary value of K[u] is equal to the eigen value 4. Therefore, the absolute 
minimum of K is the lowest eigen value Ao. 

Let the trial function be u = uo + cyuy + Colla + ..., where the {u;} are the true eigen 
functions and the {c;} are small. Then, the estimate K[u(x)] is given by 


(9.13) 


= Ao + eri + el səsə 
ltep teste 


= Age a — åo) + ely — Ao) ee (9.14) 


The results imply that K cannot be less than A, and is equal to A only if all the c; = 0, 
whence u = uo. 
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9.5.2 Example 


Estimation of the lowest eigen value 2 of 


d?u(x) 
dx? 


where u = 0 at x = 0 and x= 1. 

This equation is a Strum—Liouville equation where p(x) = 1, g(x) = 0, and p = 1. 
The analytical solution is 4 = 1°. 

Here let us estimate it with the approach described above, 1.e., the lowest eigen 
value should satisfy the following inequality: 


+ Au(x) = 0 (9.15) 


1 
/ (duldx)?dx 
— 0 


1 
/ u?dx 
0 


Consider the trial function u = 4x(1 — x), then equation gives us 2 < 10, which is 
close to x* = 9.867. Figure 9.10 shows the exact eigen function and the trial function. 

Better estimation of the lowest eigen value can be implemented by setting u 
(x)  4x(1 — x)(1 + cx) and varying the parameter, c, to minimize K. In this case, the 
‘improved’ trial function is given, its derivative is du/dx = 4(1 — 2x + 3cx* — 4cx°). 
For using EXCEL, the integrals can be approximated to 


ASK (9.16) 


1 
/ (duldx)2dx = Y\(dldx)3 Ax, 
9 i 


, È duldx); m (9.17) 
?dx = Duy Ax, and 4 — = 
f u“dx ə. x, an 


Eigen value function of u"+Au=0 
1.2 


0 0,2 0:4 0,6 0.8 1 1/2 


Trial 


Theory 


Figure 9.10. Eigen function and its trail function. 
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Table 9.4. Minimization of K[u(x)] using the trial function u(x) = 4x(1 — x)(1 + ex). 


1 c= 
2 u" u'^2 
3 i 1 1 
4 0.01 0.0099 9.8E-057 0.98 0.9604 
Rows 8 5 0.02 0.0196 0.000384” 0.96 0.9216 
6 0.03 0.0291 0.000847” 0.94 0.8836 
to 99 are 7 0.04 0.0384 0.001475" 0.92 0.8464 
r 
: 0.05 0.0475 0.002256 0.9 0.81 
hidden. isi 
100 0.97 0.0291 0.000847” -0.94 0.8836 This is the 
101 0.98 0.0196 0.000384” -0.96 0.9216 value to be 
102 0.99 0.0099 9.8E-05” -0.98 0.9604 2 
103 1 -6.7E-16 4.44E-31” -1 1 minimized. 


3.333333 N- 


a^2= 


% error = 


9.869604 
4.381084 


by taking an equal increment Ax between the integral interval [0, 1]. Table 9.4 is the 
EXCEL spreadsheet to minimize the numerical value of N/D calculated in Cell 
E105 by changing the c-value in Cell E1. The initial c-value is 0. The optimal c- 
value turns out to be 0.1325 and the eigen value, given by N/D, is 9.833, is very close 


to x = 9.867. 


Remark: in this chapter, the simplest numerical derivative and integral calculation 
are used. For their more sophisticated numerical calculations, refer to books on 


advanced numerical computation [4]. 


References 


[1] Feynman R P and Leighton R B et al 2011 The Feynman Lectures on Physics Vol. 1 (Boston 


MA: Addison-Wesley) 
[2] Goldstein 2011 Classical Mechanics (Boston, MA: Addison-Wesley) 


[3] Mathew J and Walker R L 1970 Mathematical Methods of Physics (Boston, MA: Addison- 


Wesley) 
[4] Bloch S C 2003 EXCEL for Engineers and Scientists (New York: Wiley) 


9-12 


IOP Concise Physics 


Numerical Calculation for Physics Laboratory Projects Using 
Microsoft EXCEL® 


Shinil Cho 


Chapter 10 


Quantum particle trajectories 


10.1 Objectives 


To re-formulate Nelson’s quantum stochastic dynamics for numerical computation 
and visualize trajectories of a particle in various potentials including a potential 
barrier, a harmonic potential well, and a Coulomb potential. 


EXCEL notes 


Analysis of stochastics processes requires more data points for better results. In 
order to visualize decent quantum trajectories, a set of several hundred thousand 
data points or more will be required. While EXCEL’s autofill and insert 
chart are very useful if the amount of data is relatively small, processing several 
hundred thousand data points requires an extensive amount of working time and 
labor. A quick remedy for this issue is to record the VBA macro code for the entire 
calculation and graphing process using a small data set, and then edit to expand the 
recoded macro code for more data points. Refer to appendix A.1.4 for recording 
VBA macro codes. The central limit theorem discussed in chapter 1 allows us to use 
random numbers that follow the Gaussian distribution as the quantum fluctuation. 
Appendix A.1.6 describes how to generate Gaussian random numbers using the 
NORMINV function. Menu, tool and function used in chapter 10 are: 

sv” Dropdown menu: [Insert] > [Chart] > [Scatter Chart]. 

*% Tools and Functions: Autofill, NORNINV(), RND(), COUNTIF (Where 
do you want to look?, What do you want to look for?), 
IF(Something is True, then do this, otherwise do this), 
and Macro-recording. 
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10.2 Theory —Nelson”s approach 


Nelson proposed a mathematical theory that describes the quantum mechanical 
motion of a particle from stochastic dynamics [1]. His approach is analogous to 
classical Brownian motion, and is based on Newton’s equation of motion to which 
quantum mechanical fluctuation is added. Following an article by Bacciagaluppi [2], 
we prove that Nelson’s stochastic differential equation is equivalent to Schrodinger’s 
equation. There are also research articles that prove one-dimensional cases [3, 4]. 
This book extends their work to three-dimensional systems, and draws quantum 
particle trajectories. 

Because a stochastic process is time irreversible, its forward and backward 
evolutions have a different time dependence. Nelson defines the forward and 
backward time evolutions of a quantum particle position as the following forms: 


dx() = X(t + dt) — X(t) = b,(x(t), t)dt + dii,(t) 
and (10.1) 
dx(t) = X(t — dt) — X(t) = b(x(t), dt + dw (t) 


where 5, (x(f), t) and b (x(t), t) are the forward and the backward velocities, and 
dw,(t) are Gaussian fluctuations with a mean value of 0, independent of the d¥_(s) 
for s < t (and dx_(s) for s > t), respectively. The fluctuation terms, dö, (7), follow the 
Wiener process, and their expectation values, E, at time ¢ are given by 


E| dw,(t)dw,t) | 2oöidi 
and (10.2) 
E| dw_(t)dw_(t) | = 208,dt 


where o is called the diffusion coefficient [5]. In the following discussion, we will 
show o = h/2m to derive the Schrédinger’s wave equation from the Nelson’s 
quantum stochastic dynamics. 

Because quantum mechanical trajectories are continuous and not smooth, i.e., 
non-differentiable or ‘zitterbewegung’, the velocity terms should be defined by using 
the mean forward and backward time derivatives, D, and D., as 


By(x, 1) = D,X(t) = lim —..—. X(t) = s) 

e> E 
and (10.3) 
b (x, t) = DX(t) = lim (maaa X(t) = s) 

e> E 


where (...|x(t) = x) is the expected value at time, ¢, which is conditional on the 
value. Notice that D,x(t) = D_x(t) = dx(t)/dt if X(t) is continuous and smooth as 
expected in classical physics. Figure 10.1 illustrates the forward and backward time 
derivative. 
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> time 


t-e t t+e 


Figure 10.1. Forward and backward time derivatives. 


According to Nelson, the time derivatives for an arbitrary function, f (X(t), £), 
are given by 


D.f(70), Deo = E 500 ov] re, i) 
and (10.4) 
Df (X(t), Dewar = E “57.0 = ov? se, t). 


öt 


In order to introduce the forward and backward derivatives in the Newton’s equation 
of motion, Nelson uses ‘average’ acceleration defined as @ = (D,b_ + D_b)/2. Using 
equation (10.4), the acceleration is given by 


a= əb +b.) + (E -V)b. + E -V)by - 50¥"(b, - 5). (10.5) 
Now, define 


ü = (b, - b-)/2 = (D, — D)X()/2 
and (10.6) 
7 = (5, + &)/2 = (D, + D)x(İİ2, 


and the acceleration (10.5) becomes 


a= 2a + (6 - Vö — (Ūū - V)ū — oV?i. (10.7) 

Using the acceleration (10.7), Newton’s equation of motion, mg = —V V (X), has the 
following form 

ma = Fa + (Ü - Vyö — (ü - V)u — ova = -VV(X). (10.8) 


The argument below is to prove that equation (10.8) is equivalent to the 
Sehrödinger equation (10.20). For a stochastic process, the Fokker-Plank equations 
[6] for the probability distribution, p(xX(t), t), are given by 
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2 — —V . (b,p) -- oV?p and 2 =-V- (bp) — oV?p. (10.9) 


By addition and subtraction of equation (10.9), we obtain 


1257 (bp + b_p) — Ü, i.e., 2 (3p) = 0, (10.10) 
ot 2 ot 
and 
V - (bp — bp) — oV?p =0, i.e., V- (ip — oVp) = 0. (10.11) 


According to Nelson, the solution of equation (10.11) is ip — oVp = 0, and thus 
V 
i =o" = oV(ln p). (10.12) 
p 


Taking the partial time derivative of equation (10.12) and using equation (10.10), we 
obtain 


M “uv Inp) =oV 5 =oV aly - (Gp) 
öt öt p öt p 
(10.13) 
ov(-v (0-0. Ye) 
p 
With equation (10.11), equation (10.13) becomes 

Z = ov(-v -0 -0. Ye) = —oV?0 — V(u .0). (10.14) 

p 


Define £ (77, t) = #(X, t) + (X, t) to combine the equation of motion (10.8) and 
equation (10.14) derived from the Fokker—Plank equation. Then 


E ü 1 
-i =i% €” iov? + ING - 3) + OV — BV -0 4 üV.ü — —VV 
öt at öt m 


1 
— oV?(ü + i) HV . (ü + i0)(ü + (0) — —VV 
m 


=0V E + lyg? - tyr, 
2 m 


and therefore, we obtain 


60 — (10.15) 
öt 2 m 
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Readers may notice that equation (10.15) is somewhat similar to Schrédinger’s 
equation if the gradient € term disappears. Let E = 2oV(ln T (7, t)), then the three 
&-terms in equation (10.15) become 


575 2evİ Žan n) = -2iov( 22), 
öt öt T öt 
oV?£ = 26?V?(V(in T)) = 20°¥( vər - Korp) and (10.16) 


2 


sve — vÜZeYün Tr)? - 20°9( vr] 
Therefore, equation (10.15) can be written as 


—V. rion + 202 Ly2r - Zr] =0, (10.17) 
T ot Tr m 


which means that dio 1 + lu = V depends on only time, or is only a 
t m 


function of time. 
Thus, let Dio Z + 26° vər es V = x(t), and the partial time derivative term 
m 


in this equation can be extracted to obtain 


io (X, t) m 
öt 


2 kə +++ ni t). (10.18) 
m 


t 
Now, define another function, Y(¥, t), by (x, t) = Y(X, t)exp [-2 J node |; 
With this definition, equation (10.18) becomes 
rice + 26 n = -2002 + LVF + gY. (10.19) 
öt h m 


Notice, if 2or/h = 1, i.e., o = h/2m, the y-terms in equation (10.19) are cancelled 
out. Then, equation (10.19) becomes the time-dependent Schrödinger equation: 


z 2 
Hore: t) 5 h 


— — V7 P(X, t) + VE, t). (10.20) 
ot 2m 


— t 
Because € = u + i = 2oV(InT) and I(x, t) = V(X, nexp| -2 f edil we 


obtain £ = Pia up, 
m 
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Therefore, the real and imaginary parts of £ = ï + i can be calculated by 
ü 
and (10.21) 
v= Piven Y)] = Rif S| 
m m M 


II 
| 
e) 
o 
< 
5 
=£ 
II 
sls 
e) 
o 
Fİ 1 
siz 


Back to the stochastic equation (10.1), we obtain the following formula for quantum 
particle trajectories: 


dX(t) = b,(X, t) + di, = b,(X, t)dt + je G,(t)Vdt, 

m 
and (10.22) 
dx(t) = b(X, 1) + dw. = b(X, t)dt + (2 G_(t)Vdt 


where dw, = V2o0dt G,(t) = pe G.(t)Vdt and G.(t) are normalized Gaussian 
m 


random numbers, and 


Be, ) -ü40- tfre | + wE) 
m y y 


and (10.23) 


b, t) =- +0 = tfre El + eki 
m wy wy 


The set of equations (10.23) describe the stochastic dynamics of quantum particles 
with its Schrödinger equation [6]. For simplicity, this book defines A = 1 and mass of 
particle, m = 1, in the following examples. With these units, Sehrödinger”s equation 
becomes 


W(X, t 1 
oe = - YAT, 1) + VEY, 0, (10.24) 
and the discrete step, {x;}, can be successively calculated by using 
Xil xi + Dy At + GVA (10.25) 


where the explicit forward velocity, b,, is calculated by equation (10.23) with A = 1 
and m = 1. 


10.3 Analysis of quantum particle trajectories 
10.3.1 One-dimensional free particle 


For a 1-dimensional free particle, V(x) = 0, where the wave function is a plane wave, 


i 


W(x, t) = e2". (10.26) 
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OV/0x 


Using equation (10.26), = i, and thus 


b, = hiş ə ni + In| = zl =1. (10.27) 
m wy wy 


Therefore, the discrete step, Ax; at t = ti, for the free particle is given by 
Ax; = b,,At + G(t)VAt = At + G(t)VAt, and thus 


Xin = xi + Ax; = x; + At + G(t)VAt (10.28) 


where At is the time step and G/(f) is the Gaussian random number at t = ti. 

Table 10.1 shows the initial part of the spreadsheet for calculating equation 
(10.28). In figure 10.3, there are 5,000 data points obtained by autofilling Cells 
A4: D5004. Holding the mouse to autofil1 5000 data points and drawing the 
graph is straightforward but it is time consuming. Instead, creating a VBA code to 
perform the same job is better, and EXCEL’s macro recording capability can be 
applied for this task. Figure A20 lists a macro code that is recorded for 20 data 
points, and then extended to 5000 data points after eliminating extra statements 
entered due to unintentional cursor movements. 

Figure 10.2 shows the quantum trajectory of a particle stating x = 0 at £ = 0. The 
particle is moving at the forward velocity, b, = 1, having the Gaussian fluctuation. 
Each time hitting the <F9> key, the trajectory changes due to the refreshed 
fluctuation. 


Table 10.1. Calculation of one-dimensional free particle. 


The Gaussian random number. 
“-NORM. INV (RAND (), 0,1)” 


This is the At-value. 


The initial values 
x=0 at t=0. 


The position increment: 
“=$BS1+ 
C3*SQRT ($B$1)” 


“= A3+$B$1” in -0.05 -0.24194 0.280811 0.112791; 
Cell A4. -0.12915 1.209327 0.320414 
0.191262 0.107638 0.074069 
0.2 0.265331 -1.37622 -0.25773 
8 0.25 0.007598 0.264348 0.10911 


0.116708 -0.67108 -0.10006 


Once Rovv 4 is 
completed, autofill 
Row 4 to calculate more 
data points. 


“-B3:D3” in 
Cell B4. 
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Quantum free particle 
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Figure 10.2. One-dimensionnel quantum free particle. 


10.3.2 One-dimensional potential barrier 


Suppose there is a one-dimensional potential barrier of height 1 for 0 < x < 1 in the 
above example [7]. The Schrödinger equation for this case is given by: 


OV (x,t) 1 0P (x,t) 
i = şox 


0 
öt 2 oF b 
2 
x da ain) 0<x<1 (10.29) 
ot 2 ot 
P(x) 192(x,0) 
= $ ie 
yr 2 or . 


By setting P(x, t) = y/((c)e” Fİ, the space part, y(x), is given by 
"A69k ə x <0 


+i —ikx ky =V2mE 
vv, (9) =C +D; O<x<1 where | TY" and V, =1. (10.30) 
k = ,i2m(E -V,) 


tiko . 


y,(x)=T-e"; xol 


Table 10.2 lists the forward velocity terms when £ = 0.5 (1.e., E < Vo). The 
spreadsheet for this example is similar to the one shown in table 10.2. However, 
because the wave function is different for a different coordinate range, the condi- 
tional position increment, Ax; at t = ti, must be selected depending on where the 
particle is. Notice that EXCEL has an IF(Something is True, then do 
this, otherwise do this) function, which can be used to select appropriate 
velocity term, depending on the current particle position. 

Table 10.3 is the initial part of the spreadsheet. In this case, we took a short time 
increment (0.0001 instead of 0.005) and 50 000 data points. After completing the 
Ax-calculation, Columns E and F are added to obtain the position distribution using 
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Table 10.2. Forward velocity for a potential barrier. 


x bi(x, t) 
x<0 1 — (tanh(1))? + 2(tanh(1))cos(2x) 
1 + (tanh(1))sin(2x) 
— tanh(2x) + ——— 
cosh(2x) 
x> 1 1 


Table 10.3. Spreadsheet configuration for the potential barrier problem. 


The position increment Ax is given by 
= (IF (B3<0, ((1-$D$1*2+2*$D$1*COS (2*B3) ) / (14+$D$1*2+2*$D$1*SIN (2*B3) )) , 0) 
+IF (AND (B3>=0, B3<=1) , TANH (2*B3) + (1/COSH (2*B3) ) , 0) 

+IF (B3>1,1,0)) *$B$1+C3*SQRT ($B$1) 


Cell D1 stores tanh(1). 


F c| Taking the 


1 0.0001 -0.76159 tally for each 
2 time x A(t) Ax(t) we 

3 0 -2 -0.86976105 -0.008366301 -10 0 position 

4 0.0001 -2.008366301 0.623002862 -0.008357724 -9.995 0 interval using 
5 0.0002  -2.016724025 -0.102638685 — 0.006578841 -9.99 0 


COUNT IF. 


the COUNT IF function described in section 1.4.2. Let the possible range of the 
position taken be [10, +10] with the step 0.005 in Cells E3 : E4003. The frequency 
(tally) counting statements are similar to what we used in section 1.4.2. For this 
example, they are 

=COUNTIF ($B$3:$B$50003,"="&E3) in Cell F3, and 

=COUNTIF ($B$3:$B$50003,"<="&E4) -COUNTIF ($B$3:$B 
$50003,"<="&E3) in Cell F4, and so forth. 

Figure 10.3 shows a typical trajectory where the region of the potential barrier, 0 
€ x < l, is highlighted and figure 10.4 is the position distribution of the particle 
across the potential barrier. Tunneling due to the quantum fluctuation is clearly 
noticeable. 
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One-dimensional potential barrier 


Figure 10.3. Particle trajectory tunneling a potential barrier. 
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Figure 10.4. Position distribution across the potential barrier. 


10.3.3 Harmonic oscillator 
The Schrödinger equation for a particle in a one-dimensional harmonic potential is 
given by 

jae = .. t)+ .... t) (10.31) 


where the spring constant is k = 1. Table 10.4 lists the wave functions and the 
forward velocity terms for several lower energy levels. 
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Table 10.4. Wave functions and forward velocities of the harmonic oscillator. 


Energy level: 
En = n+1/2, 
nz0,1,2,3, ... Wave function bx, t) =b 
1 UE 121 -—x 
Eo = 2 P(x, t) = (2) enə 
T 
3 1211/4 Lə 1) 
B= = EE D ——x”—i—t b--İx —— 
1 7 P(x, t) (5) 2) 2xe 27772 R 
5 1271 \"4 Lo 6s ( 4x ) 
E=2> x palt) (2) ax- per poets 
2 P(x, t) (3) (+) (4x4 — 2)e 27 72 me] 
7 1 1/2 1 1/4 122 47 6x2 — 3 
E= 5 va, o = (5 (2) 8x3 — 12x)e72* “iə b = -|x - — 
2 Ye) 48 m - me ` 2x3 — 3x 
9 i 1x22 8x(2x2 — 3) 
E= > W(x, t) = (Gi) (=) 16x4 — 48x? + 12)e27””7f b = -| x — — —”. 
2 227068 ) 4x4 — 12x? + 3 


The spreadsheet for the harmonic oscillator is similar to the one shown in 
table 10.3. Only the expression of the position increment, Ax, needs to be changed. 
The COUNT IF function is also used to find the probability distribution of the 
position in the position range [—6, +6]. Figure A21 lists the VBA macro code 
recorded for the ground state (n = 0 and 50 000 terms) including the position 
probability calculation. Notice that the VBA code was created in the follow manner: 
(i) a macro code was recoded after completing the trajectory part (with 100 terms); 
(ii) the code was merged with the trajectory part; and (iii) the number of data points 
of the merged code were changed to 50 000. 

Figure 10.5 shows the trajectories, and the position distributions of the harmonic 
oscillator. The point distribution simulates the probability distribution from the 
wave function. In these analyses, the numbers of data points increase from 50 000 
for the n = 0 state to 200 000 for the n = 4 state. The time increment is 0.01 for all 
cases. The occasional spikes are noticeable when n is larger. They may occur when 
the forward velocity values take large values. The time interval and the number of 
data points are determined empirically. 
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Figure 10.5. Trajectories of a harmonic oscillator. 
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2500 


Figure 10.5. (Continued.) 


10.3.4 Hydrogen atom 


The Schrédinger equation for a hydrogen atom is given by 
1 = ... t) + lyg, t) where r=Ę4x? +y? +z? (10.32) 
r 


öt 
where we set e7/4ze) = 1 for the Coulomb potential. Table 10.5 lists the wave 
functions of the lower energy levels and the forward velocity terms. 

Since the electron orbit in a hydrogen atom is three-dimensional, each of the 
three-dimensional coordinates has an independent Gaussian random number, i.e., 
there are three independent sets of Gaussian random numbers, A(t), B(t), and 
C(t), as shown in table 10.6. 

Figure 10.6 shows the two-dimensional projections of the electron trajectories of a 
hydrogen atom. These graphs are selected from several trajectories after running the 
macro codes several times by hitting the <F9> key. For n = 1 and 2, we used the xy- 
plane and for n = 3, we used the xz-plane. Notice the coordinate scales of the 
trajectories are different for different energy states because the ‘drifting’ region of the 
electron becomes larger as the principle quantum number increases. The number of 
data points used for these trajectories are 200 000 for the n = 1 state, 600 000 for 
n=2 states, and 1 000 000 for the n = 3 states! Figure A22 lists the VAB code for the 
n=3, £ = 2, and m, = 0 state. 


Remark: readers should remember that the computational performance depends on 
the hardware configuration of the computer used. For example, for these calcu- 
lations, the author used a home-built PC that has a quad CPU (3.20 GHz), 16GB of 
RAM, a SSD, and a graphics card with 2GB RAM. A notebook PC with a quad 
CPU of 3.5 GHz with 83GB RAM was unable to complete the 1 million data point 
calculations due to a lack of memory. It is really amazing that a mediocre home- 
built PC can handle the computations of a vast number of data points while also 
accessing the Internet! 
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Table 10.6. Initial part of the spreadsheet for hydrogen atom (n=1). 


R(t) calculates r=SQRT(x*2+y*2+z42), 


There are 3 independent sets of Gaussian 
random numbers for x, y, and z at each time. 


Ax is given by 
=- (E3/$H3) *$B$1+B3*SQRT 


B(t) C(t) Y(t) Z(t) R(t) BeltaY Delta Z 

o! -1.30525 1.51754 -0.469061 1 1 1 1.732051 -0.136 0.145981 -0.05268 
0.01! -0.83448 -0.6882 -1.70952 1 0.863702 1.145981 0.94732 1.719496 -0.08847 -0.07548 -0.17646 
0.02! -0.20482 ..0.13923 _-1.70426/ 0.775231 1.070496 0.770859 1.530087 -0.02555 .0.02097 -0.17546 


H-atom (100) 


6D 8 


-6 


Figure 10.6. Orbital electron trajectories in a hydrogen atom. 
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Figure 10.6. (Continued.) 
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H-atom (211) 
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Figure 10.6. (Continued.) 
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Figure 10.6. (Continued.) 
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Figure 10.6. (Continued.) 
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Figure 10.6. (Continued.) 
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Figure 10.6. (Continued.) 
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Appendix A 


A.1 EXCEL options 
A.1.1 Autofill 


Autofill is a very useful EXCEL tool that is frequently used for scientific 
calculations. With its capability, iterative calculations can be carried out without 
programming. Here is a simple example that enters integer 0, 1, 2, 3, ... in 
Column A: 

(1) Input 0 in Cell A1. 

(2) Enter =A1+1 in Cell A2 and press <Enter>. The value of Cell A2 
becomes 1. 

(3) Place the cursor in Cell A2 and hit <Enter>. The cell is emphasized. 

(4) Click on the fill handle (a small solid square at the lower right corner). The 
cursor becomes a plus sign (+). While pressing the right mouse button, drag 
the cursor to the cells below in Column A until reaching the desired integer 
value (figure A1). 


| Enter 0 | 


o ] Enter =A1+1 | 


1 
2 
3 
4 
5 
6 
m 
8 


1 
2 
3 
4 
5 
6 
7 
8 


Drag [+] toa 
desired row 


Figure Al. Autofill. 
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A.1.2 Adding ‘data analysis’ 


Because EXCEL’s default setting is not made for scientific computations, we need to 
add [Data Analysis Tools] in order to add the histogram option and other 
useful scientific calculation tools [1]. From the [File] menu, select [options] to 
display the “EXCEL Options’ screen. Click on [Add-ins] to display the following 
screen (figure A2). 

Next, click on [Go] to display available add-ins, and then check [Analysis 
ToolPakl [Analysis ToolPak - VBAİ, and [Solver Add-in] (figure A3). 


A.1.3 Enabling VBA macro 


EXCEL macro is a Visual Basic (VB) programing environment. Appendix A3.4 lists 
the source codes created for this book. Those who are interested in EXCEL macro, 
refer to ‘Excel? VBA for Physicists’ of this Concise Series [2]. Take the following 
steps to enable EXCEL’s macro capability: 
(5) Go to [Trust Center]. 
(6) From [Option], go to [Trust Center], and click on [Trust Center 
Settings..] of [Excel Options] in [Microsoft Excel Trust 
Center] (figure A4). 


Excel Options 


ced View and manage Microsoft Office Add-ins. 


Location 


Active Application Add-ins 
Analysis ToolPak 


Inactive Application Add-ins 
Address (Korean Address Recognizer 


Figure A2. EXCEL options Screen. 
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Add-ins 


Add-ins available: 


Hİ Analysis ToolPak - VBA 


LT Euro Currency Tools Cancel 
Kİ Solver Add-in 


Brovvse... 


Automation... 


Solver Add-in 


Tool for optimization and equation solving 


Figure A3. Available Add-ins. 


Excel Options 


G I 
ə © Help keep your documents safe and your computer secure and healthy. 


Formulas 


Proofing 


Security & more 


Save Visit Office.com to learn more about protecting your privacy and security. 


Language Microsoft Trustworth mputin 


Ease of Access 
Microsoft Excel Trust Center 
Advanced 


The Trust Center contains security and privacy settings. These settings help keep your computer 
Customize Ribbon secure. We recommend that you do not change these settings. Trust Center Settings... 
Quick Access Toolbar 
Add-ins 


Trust Center 


Figure A4. EXCEL options for enabling macro. 
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(7) Select [Macro Settings] and check [Enable all macro (not 
recommended; potentially dangerous code can run) ]in order 
to use this capability. Click [OK] to complete the setting (figure A5). 

(8) In the pulldown menu, click on [View] and click on [Macros] to select 
[View Macros] and create a macro program (figure A6). 

(9) After entering a macro name, we can create its source code using a built-in 
editor (figure A7). 


Trust Center 


Trusted Publishers R 
Macro Settings 
Trusted Locations 
O Disable all macros without notification 
Trusted Documents 
ith notification 


ed Add-in Catalogs 
3 macros 


Add-ins @ Enable all macros (not recommended; potentially dangerous code can run) 


ActiveX Settings ? 
Developer Macro Settings 
Macro Settings 


Trust access to the VBA project object mode 


Protected View 
Message Bar 
External Content 
File Block Settings 


Privacy Options 


Figure A5. Macro settings. 


Macro ? x 


Macro name: 


il ” Run 
Step into 
_\ Edit 
Enter your macro name Create 
and click on [Create]. | Delete 
l Options... 
Macros in: | All Open Workbooks {v 


Description 


Cancel 


Figure A6. Macro entry. 
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# Microsoft Visual Basic for Applications - Book! - (Module! (Code)) - o x 

A$ Ele Edit View Insert Format Debug Bun Tools Add-ins Window Help - x 
@@-a A” » n a M By YY ? O Waco BE 

Project “VBAPTOyeet Xİ İ Gem Je ox 
u iz ş Sub LP() “al 


E Bİ atpvbəenixis (A1 ^ End Sub — 


Write your source code here. 


Sheeti (St 
3 Thisworkb 
5 (fiy Modules 
AR Modulel 
© BÍ veaproject (Fun Y 
< > 


Figure A7. Macro editor. 


© Formula Bar 


Normal Page Break Page Custom | (7) gridlines [Z Headings Zoom 100% Zoomto Al 
revi s 


range Fi 
lev: Layout View: Selection All Panes ~ 


Workbook Views Show Zoom 


AL x fe 


Figure A8. Preparing a macro recording. 


A.1.4 Recording macro code 


Recording macro is a quick way to make a macro code, which can be saved and 
edited for other data analysis. The following example is to ‘draw a scatter graph of 
y = sin(x) + 2cos(x) for x = —z to + z by 0.1 step.’ Refer to section A.1.3 for enabling 
macro functions: 

(1) Open a new sheet. 

(2) Click on the [View] pull down menu and then click on the downward 
triangle (W) of [Macro] to select [Record Macro...] (figure A8). 

(3) Enter a macro name and hit [OK] (figure A9). 

(4) Make sure that [Stop Macro] will be displayed instead of [Record 
macro... in the [Macro] menu but do not click on it yet. Close the 
[Macro] menu by moving the cursor elsewhere and click there (figure A10). 

(5) Enter x in Cell A1, and y in Cell B1. Enter =-pi() in cell A2, and =A2 
+0.1 in Cell A3. Autofill Cell A2 to Cell A70 where the x-value exceeds 
+n. Delete Cells A66 to A80. Enter =sin(a2) *2*cos(2*a2) in Cell B2. 
Autofill Cells B2:B65. 

(6) Highlight Cells A2 :B65 and insert [Scatter with smooth line with 
markers]. Figure All shows the graph of data of Row 2 to Row 62 
without the chart and axes titles. 

(7) Click on the [Macro] menu and click on [Stop macro] to stop recording 
the macro. Click on the [View Macro] of the [Macro] menu to run or edit 
the created macro code. Select the [Edit] to see the created macro code 
(figure A12). 
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Record Macro 


Macro name: 


graphingTrigfunction 


_ Shortcut key: 
Ctrl+ 


Store macro in: 


This Workbook 


Description: 


Insert Page Layout Formulas Data Review 


(78 Bo E Formula Bar Q Lə cS EE isp msl 
əzə, = = ream! zə Hide 

Normal Page Break Page Custom (7) gridlines [Z] Headings Zoom 100% Zoom to | New Arrange Freeze Switch 

Preview Layout Views Selection Window All Panes ~ Windows ~ 


Workbook Views Show Zoom 


Tse Relative References 


Figure A10. Stopping recording macro code. 


Chart Title 
3 


Figure A11. Graph created while recording macro. 
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Macro 


Macro name: 


graphingTrigfunction 


Step Into 


Create 
Delete 


Options... 


Macros in: All Open Workbooks 


Description 


Figure A12. Macro code menu. 


(8) Figure A 13 lists the created macro code. Reflecting what was entered in the 
sheet, it is possible to roughly think about what the codes do. 


A.1.5 Enabling iterative calculation 


From the [File] menu, select [options] to display the [EXCEL Options] screen. 
Click on [Formulas] to display the following screen (figure A14). 


A.1.6 Generating Gaussian random numbers 


EXCEL does not have a built-in function that generates Gaussian-distributed 
random numbers. However, combining two EXCEL functions can generate random 
numbers that follow the Gaussian distribution of the mean = 0 and the standard 
deviation = 1 (‘Gaussian random number’): 
(a) RAND() returns a uniformly distributed random number that is greater 
than or equal to 0 and less than 1; and 
(b) NORM. INV () returns the inverse of the normal cumulative distribution for 
the specified mean and standard deviation. 


Thus, combining these functions together, NORM.INV(RAND(),0,1) generates 
Gaussian random numbers. 
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Sub grap! 


End Sub 


" GraphingTrigfunction Macro 


hingTrigfunction() 


Application. WindowState = xIMinimized 
Application. WindowState = xINormal 
ActiveCell.FormulaR1C1 = "x" 
Range("B1").Select 

ActiveCell.FormulaR1C1 = "y" 
Range("A2").Select 

ActiveCell.FormulaR1C1 = "—-PI()" 
Range("A3").Select 

ActiveCell.FormulaR1C1 = "=R[-1]C+0.1" 


Range("A3").Select 

Selection. AutoFill Destination:=Range("A3:A80"), Type:=xlFillDefault 
Range("A3:A80").Select 

ActiveWindow.SmallScroll Down:=-63 


These 3 lines describes the 
action of deleting extra 
Cells A66 to A80. Notice 
these lines can be 
eliminated if A80 is 
changed to A65 in the 
stamen above. 


Selection.ClearContents 


Range("A66:A80").Select 
ActiveWindow.SmallScroll Down:=-69 + 


Range("B2").Select 

ActiveCell.FormulaR1C1 = "=SIN(RC[-1])+2*COS(2*RC[-1])" 
Range("B2").Select 

Selection. AutoFill Destination:=Range("B2:B65"), Type:=xIFillDefault 
Range("B2:B65").Select 

ActiveCell.FormulaR1C1 = "=SIN(RC[-1])+2*COS(2*RC[-1])" 


Range("A2:B65").Select 
ActiveSheet.Shapes.AddChart2(240, xIXYScatterSmooth).Select 
ActiveChart.SetSourceData Source:=Range("Sheet1!$A$2:$B$65") 
ActiveSheet.Shapes("Chart 1").IncrementLeft -169.1320472441 
ActiveSheet.Shapes("Chart 1").IncrementTop -16.3018897638 


Range("R54").Select 
ActiveSheet.ChartObjects("Chart 1").Activate 
ActiveSheet.ChartObjects("Chart 1").Activate 
ActiveChart.ChartArea.Copy 
Application. WindowState = xIMinimized 
Application. WindowState = xINormal 
Application. WindowState = xIMinimized 
Application. WindowState = xINormal 

Range("R50").Select 


Figure A13. Recorded VBA macro code. 


Figure A15 shows an example of distribution of random numbers generated by 
NORM. INV (RAND () ,0,1), and the tally of frequencies in an interval [—60, +60] 
where <x> = 0 and o = 1. 


(1) 
(2) 
(3) 


Enter = NORM.INV(RAND() ,0,1) in Cell A2. Use autofil1 to generate 
1000 Gaussian random numbers to Cell A1001. 

Enter =-6 in Cell C2, and =C2+0.5 in Cell C3. Use autofill to generate 
-6 to +6 by step 0.5 in Column C. 

For taking a tally in Column D, we used the COUNTIF function used in 
chapter 1. The COUNTIF function counts the number of cells that meet a 
criterion, e.g., to count the number of times a particular value appears in 
Column A. 

In this case, enter =COUNTIF ($A$2:$A$1001, "<="&C3) in Cell D2 
to find how many random numbers in Column A are less than and equal to 
the value of Cell C3, 1.e., the number of occurrences of the Column A that 
are <6, and enter =COUNTIF($A$2:$A$1001, "<="&C4) -COUNTIF 
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Excel Options ? x 


General ə) n ı 
Fc-. ki fx Change options related to formula calculation, performance, and error handling. 

Proofing Calculation options 

Save Workbook Calculation G O Enable iterative calculation 

Language © Automatic Maximum Iterations: | 100 2 


Automatic except for data tab) 
Ease of Access á - Maximum Change: [0.001 


Advanceg 


Customi Check this box to enable iterative 
Quick Aq calculation. Also determined 
Add-ins Maximum iterations, e.g., 


10000 and Maximum change, 
e.g., 0.000001. 


m Reset Ignored Errors 
Indicate errors using this color: | £2 7 


Error checking rules 


Cells containing formulas that result in an error & Formulas which omit cells in a region © 
İZ) inconsistent calculated column formula in tables < Unlocked cells containing formulas © 
İZ) Cells containing years represented as 2 digits € C Formulas referring to empty cells © 

El Numbers formatted as text or preceded by an apostrophe © Data entered in a table is inyalid ‘ 

İZİ Formulas inconsistent with other formulas in the region < Misleading number formats © 


Figure A14. Enable iterative calculation. 


A 

Norm.INV 
0.272594675 y (frequency) 
1.267308917 
-0.056375743 
0.563126195 
-0.008612423 
0.898059444 
-1.140388283 
-0.56031151 
2.00749483 
0.563956656 n 
-0.808288564 i m 
-0.056807813 
-1.701376021 
-0.291208286 
0.814246343 
-0.049902004 
-0.427645644 


gogoxzococscoo 


. . 
.........22:....2728 


6 -5 4 -3 -2 -1 0 1 2 3 4 5 6 


Figure A15. Gaussian random number. 
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($A$2:$A$1001,"<="&C3) in Cell D3 to find how many random 
numbers in Column A are less than and equal to the value of Cell C4 and 
how many random numbers in Column A are less than and equal to the 
value of Cell C3, i.e., the number of occurrences of the Column A that are 
5.5 < [number] < 6. Autofill to generate the criterion to Cell D1203. 


A.2 Calculation of firing speed of a rolling ball 


The gravitational potential energy of a ball at rest at the starting point from the 
firing position, APE, is equal to the kinetic energy of the ball, AKE, at firing point: 
APE = AKE. Because the ball is rolling down (without slipping) on the guiding rail, 
the kinetic energy originates from the translational and the rotational motions, 


1 
AKE = [AKE translation] + [AKE otation] where AKErranslation — —muq. 


2 
The rotational kinetic energy should be given by 
2 
AKE rotation = las = “ED = Lae 
2 2\5 5 r2 
— —muğ e Liu l - 
05:7R)—2 5 de 


where R is the radius of the ball, 2w is the rail width, m is the mass of the ball, J = 
(2/5)mR”, the angular velocity, œ = vp/r as shown in figure A 16. Thus, we obtain the 
total kinetic energy 


1 1 
AKE = —muş + 50 —  /uyı” (A.2) 
L- (x) 


For our equipment, w/R = 0.635, and thus 


vo = Jgh. (A.3) 


A rolling ball 


a ə 


.— 
|] A rail track 


SSS 


Figure A16. A ball rolling fired at the projectile launcher. 
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A.3 RLC oscillator circuit 


Figure Al7 shows an RLC circuit containing a nonlinear amplifier with positive 
feedback. Suppose the nonlinear amplifier outputs current, i,, has a form f(V), where 
V is the input voltage, which is a small signal voltage, v,(7), on the bias voltage, Va: 
io =f (Vg + 2). 

Let the bias voltage Ve set to the operating point where the small signal gain takes 
the maximum value, then to the first nonlinear term, 


; of V) 9707) 3 3 
lo = Vg + v(t = < v(t) + Us = 80s — v(t A.4 
Petan = y m Ut- rA .. gü, — 8305(t) (A.A) 
2 3 
where g, = FSV) ə. - or), = 0, and g, = 2) > 0 at the 
ÖV İy-v, OV? Wer, öV? Hn, 


operating point. 
The voltage across the circuit, v, is given by 


dir(t) 


ot) -L dı 


= Rig(t) = = ic(t)dt and i(t) + ig(t) + ic(t) + i,(t) = 0. (A.5) 


Therefore, 


d z ; : _ v(t) 1 do(t) dolt) didt) _ 
OO = RO Pic) 4) = R. ui FT +C a + Ai =0. (A.6) 


Suppose v, (t) is a positive feedback, v, = —fv, then ip = -gv + gv? from 
equation (A.4) and equation (A.6) becomes 


v(t) 1 do(t) dlt) ( do 2) 
— + — +C + | -g P—  3g: 9502 | 0, A.7 
L Rat dh, urn." r” -. 
and therefore, 
2 2 
awd _ ly ra (A.8) 
dit, dt LC 


Figure A17. RLC oscillator circuit. 
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where 
R-1 3g:93 
a= SPE E 8P l 
RC C 
Let v = ECM then 
£ 
d?x dx 1 
— — e(l — x? + — 0 
dt? ... Vit ıc” 
Finally, let t = V LCz, and we obtain the van der Pool equation: 
d?x dx 
— “gü — x*)— + x =0 
uu: 


where 


w= aI = |E(ep- x) 


A.4 Another circuit for exercising Kirchoff’s law 


(A.9) 


(A.10) 


(A.11) 


(A.12) 


Figure A.18 shows a circuit where Ry = 1.0 KQ, R3 — 2.2 kO, R, = 3.8 kO, Rs = 2.2 kO, 
Re = 1.5 kQ, and Vo = 6.3 V. Let us calculate the currents in mA across each resistor. 


From the Kirchhoff”s Law, we obtain 


-bh-h=0 at the junction a, 
h-h-k=0 at the junction 5, 
h+h-k=0 at the junction e, 

Roh + Rolo = V for the loop abcdf, 
Rab + Rsls = V for the loop aecdf, and 
Rab + Ra t+ Ri = V for the loop abecdf 


Figure A18. An example of a circuit. 
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where we did not include 7, = I5 + J¢ because it is not an independent equation. 
Using the resistors and the battery values, the set of above equations become the 
matrix form, AI = V: 


1-1-1 0 0 0111 0 
10 0 0 -1 -1] İR 0 
o1 0 -1-1 0) 16). |o 
01 0 022 0} İLİ 163) 
00 22 0 0 15] |Z| 163 
01 0 38 0 15] [4s] 163 


The current vector can be calculated by 7 = A~'V. The reader should try to find the 
final result of figure A19. 


10 | 0.543031 0.456969 0.190784 0.294433 0.246832 0.057813 
11| -0.12567 0.125666 0.464966 0.268469 -0.05712 0.140899 


12| -0.3313 0.331302 -0.27418 0.025964 0.303953 -0.08309 
13| -0.18279 0.182788 -0.32369 -0.06404 -0.08309 0.204944 
14| 0.057121 -0.05712 -0.21135 0.332514 0.025964 -0.06404 

0.48591 -0.48591 0.402133 -0.03808 0.220868 0.121858 


7.548397 
4.438309 
3.110088 

0.72845 
3.709859 
3.838538 


Figure A19. Result of 7 = A7'V. 


A.5 Macro codes for quantum particles 


One-dimensional free particle 

Figure A20 is the macro code for 5000 data points of a one-dimensional free 
particle by modifying a recoded macro code for 20 data points. In this program, only 
the numerical values ‘5000’ in red were changed from ‘20.’ 

One-dimensional harmonic oscillator in the ground state after changing the number 
of data points from ‘100’ to ‘50 000” (figure A21). 

Hydrogen atom in the (n = 3, 1 — 2, m, = 0) state after changing the number of data 
points to 1000 000 (figure A22). 
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Range ("A3") .Select 
ActiveCell.FormulaR1C1 = 
Range ("B3") .Select 
ActiveCell.FormulaR1C1 = 
Range ("C3") .Select 
ActiveCell.FormulaR1Cl = 
Range ("D3") .Select 
ActiveCell.FormulaR1C1 = 
Range ("D3") .Select 
ActiveCell.FormulaR1C1 = 
Range ("A4") .Select 
ActiveCel1.FormulaR1C1 = 
Range ("B4") .Select 
ActiveCell.FormulaR1C1 = 
Range ("C4") .Select 
ActiveCell.FormulaR1C1 = 
Range ("D4") .Select 
ActiveCell.FormulaR1Cl = 


Range ("A4:D4") .Select 


Selection.AutoFill Destination:=Range ("A4:D5000"), 


Range ("A3:B5000") .Select 


ActiveSheet.Shapes.AddChart2 (240, 


non 
non 
"-NORM. INV (RAND () , 0,1)" 
"$B$1+SQRT ($B$1) *C3" 
"=R1C2+SQRT (R1C2) *RC[-1]" 
"=SR[-1]C+R1C2" 
"=R[-1]C+R[-1]C[2]" 
"-NORM. INV (RAND() ,0,1)" 


"-RİIC2:3ORT (R1C2) *RC[-1]" 


Type :=x1lFillDefault 


xlXYScatterLinesNoMarkers) .Select 


ActiveChart.SetSourceData Source:=Range ("Sheet1!$A$3:$B$5000") 


End Sub 


Sub Freeelectron() 
"screen deletion Macro 
Cells.Select 
Selection.Delete Shift:-xlUp 
"Quantum free particle 
Range("Al").Select 
ActiveCell.FormulaR1C1 = 
Range ("B1") .Select 
ActiveCell.FormulaR1C1 = 
Range ("A2") .Select 
ActiveCell.FormulaR1Cl = 
Range ("B2") .Select 
ActiveCell.FormulaR1Cl = 
Range ("C2") .Select 
ActiveCell.FormulaR1C1 = 
Range ("D2") .Select 
ActiveCell.FormulaR1C1 = 


"Delta t=" 
POSS 
"Time" 

"x" 


"A(t)" 


"Delta(t)" 


Figure A20. Macro code for a one-dimensional quantum free particle. 
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Sub Harmonic0() 
" screendeletion Macro 
Cells.Select 
Selection. Delete Shift:=xlUp 
" Quantum harmonic oscillator (n=0) 
Range ("Al") .Select 
ActiveCell.FormulaR1C1 
Range ("B1") .Select 
ActiveCell.FormulaR1C1 = "0.005" 
Range ("A2") .Select 
ActiveCell.FormulaR1C1 = "Time" 
Range ("B2") .Select 
ActiveCell.FormulaR1C1 = "X(t)" 
Range ("C2") .Select 
ActiveCell.FormulaR1C1 = "A(t)" 
Range ("D2") .Select 
ActiveCel1.FormulaR1C1 = "Delta X" 


1565. 


Range("A3").Select 
ActiveCell.FormulaR1C1 = "0" 
Range ("B3") .Select 
ActiveCell FormulaRlel = "0" 
Range ("C3") .Select 
ActiveCell.FormulaR1C1 = "=NORM.INV(RAND(),0,1)" 
Range ("D3") .Select 
ActiveCel1.FormulaR1C1 = "=—-RC[-2]*R1C2+SQRT(R1C2)*RC[-1]" 


Range ("A4") .Select 

ActiveCell.FormulaR1C1 = "=R[-1]C+R1C2" 
Range ("B4") .Select 

ActiveCell.FormulaR1C1 = "=R[-1]C+R[-1]C[2]" 
Range ("C3") .Select 

ActiveCell.FormulaR1C1 = "=NORM.INV(RAND(),0,1)" 


Range ("C3:D3") .Select 
Selection.AutoFill Destination:=Range("C3:D4"), Type:-xlFillDefault 


Range("A4:D4").Select 
Selection.AutoFill Destination:“Range("A4:D50000"), Type:=xlFillDefault 


Range ("A3:B50000") .Select 
ActiveSheet.Shapes.AddChart2 (240, x1XYScatterLinesNoMarkers) .Select 
ActiveChart.SetSourceData Source:=Range ("Sheet1!$A$3:$B$50000") 


“ Calculation of position probability 
Range ("E2") .Select 


ActiveCell.FormulaR1C1 € "Bin" 
Range("F2").Select 

ActiveCell.FormulaR1C1 — "Freq" 
Range("E3").Select 

ActiveCell.FormulaR1C1 = "-6" 


Range("F3").Select 
ActiveCell.FormulaR1C1 

Range ("E4") .Select 
ActiveCell.FormulaR1C1 

Range("F4").Select 
ActiveCell.FormulaR1C1 “r 
"-COUNTIF (R3C2:R50000C2, ""c-""SRCİ-11)-COUNTİF (R3C2:R50000C2," 


"=COUNTIF (R3C2:R5000062, ""xe"",RCİ-11)" 


uzalı 


"SRI-11C1-11)" 


Range("E4:F4").Select 
Selection.AutoFill Destination:”Range("E4:F1203"), Type:-xlFillDefault 


Range("E3:F1203").Select 
ActiveSheet.Shapes.AddChart2(240, xlXYScatterLinesNoMarkers) .Select 
ActiveChart.SetSourceData Source:-Range("Sheet11$E$3:$F$1203") 
End Sub 


Figure A21. VBA code for the ground state harmonic oscillator (50 000 terms). 
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Sub H320xzy() 


screendeletion Macro 
Cells Select 
Selection.Delete Shift:=xlUp 
Hydrogen atom (320-state) 
Range ("Al") .Select 
ActiveCell.FormulaR1€1 = "Delta t =" 
Range ("B1") .Select 
ActiveCell.FormulaR1C1 = "0.02" 
Range ("A2") .Select 
ActiveCell.FormulaR1Cl = "t" 
Range ("B2") .Select 
ActiveCell.FormulaR1C1 = "A" 
Range ("C2") .Select 
Activecell- Frormulari¢l = "B" 
Range ("D2") .Select 
ActiveCell.FormulaR1C1 = "C" 
Range("E2").Select 
ActiveCell.FormulaR1C1 = "X(t)" 
Range("F2").Select 
ActiveCell.FormulaR1C1 = "Z(t)" 
Range ("G2") .Select 
ActiveCell.FormulaR1C1 = "Y(t)" 
Range("H2").Select 
ActiveCell.FormulaR1C1 = "R(t)" 
Range("I2").Select 
ActiveCell.FormulaR1C1 = "D(t)" 
Range("U2").Select 
Activečcell. FormulaRicr = "Delta xK" 
Range ("K2") .Select 
ActiveCell.FormulaR1C1 = "Delta Z" 
Range ("L2") .Select 
ActiveCel1.FormulaR1C1 = "Delta Y" 


Range ("A3") .Select 
ActiveCell.FormulaR1C1 = "0" 
Range ("B3") .Select 
ActiveCell.FormulaR1C1 = "=NORM.INV(RAND(),0,1)" 
Range ("C3") .Select 
ActiveCell.FormulaR1C1 = "=NORM.INV(RAND(),0,1)" 
Range ("D3") .Select 
ActiveCell.FormulaR1C1 = "=NORM.INV(RAND(),0,1)" 
Range ("E3") .Select 
ActiveCell.FormulaR1C1 = "0.1" 
Range ("F3") .Select 
ActiveCell.FormulaR1C1 = "1" 
Range ("G3") .Select 
Actiyecell, FormulaRici = Ers 
Range ("H3") .Select 
ActiveCell.FormulaR1C1 = "=SQRT(RC[-3]%*2+RC[-2]*2+RC[-1]%*2)" 
Range ("I3") .Select 
ActiveCell.FormulaR1C1 = "= 3*RC[-3]*2-RC[-1]*2" 
Range ("J3") .Select 
ActiveCell.FormulaR1Cl = _ 
u—  (BGTE51/(2“RO 1221 )5:2 RO: /RC[-1]] *RICZ4RO[-7] "SORT (5162) 
Range("K3").Select 
ActiveCell.FormulaR1C1 = _ 
n= (RG ON (25RCT 3 10): RG EE İ//EG L272) RIC24RC LEO5SOET (RICZ) ee 
Range("L3").Select 
ActiveCell.FormulaR1Cl = _ 
"= (RC[-5]/ (2*RC[-4])+2*RC[-5]/RC[-3])*RIC2+RC[-8]*SORT(RIC2)" 


Range ("A4") .Select 


ActiveCell.FormulaR1C1 = "=R[-1]C+R1C2" 
Range ("B4") .Select 

ActiveCell.FormulaR1C1 = "=NORM.INV(RAND(),0,1)" 
Range ("C4") select 

ActiveCell.FormulaR1C1 = "=NORM.INV(RAND(),0,1)" 
Range ("D4") .Select 

ActiveCell.FormulaR1C1 = "=NORM.INV(RAND(),0,1)" 
Range ("E4") .Select 

ActiveCell.FormulaR1C1 = "=R[-1]C+R[-1]C[5]" 
Range ("F4") .Select 

ActiveCell.FormulaR1C1 = "=R[-1]C+R[-1]C[5]" 
Range ("G4") .Select 

ActiveCell.FormulaR1C1 = "=R[-1]C+R[-1]C[5]" 


Range ("H3:L3") .Select 
Selection.AutoFill Destination:=Range("H3:L4"), Type:-xlFillDefault 


Figure A22. Hydrogen atom 320-state. 
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Range ("A4:L4") .Select 
Selection.AutoFill Destination:“Range("A4:L1000000"), Type:exlFillDefault 


Range ("E3:F1000000") .Select 
ActiveSheet.Shapes.AddChart2 (240, x1lXYScatterSmoothNoMarkers) .Select 
ActiveChart.SetSourceData Source:=Range ("Sheet1!$E$3:$F$1000000") 


ActiveChart.ChartTitle.Select 
Selection.Caption = "H-atom (320)" 


ActiveChart.FullSeriesCollection (1) .Select 

With Selection.Format.Line 
-Visible = msoTrue 
-Weight = 0.25 
.ForeColor.ObiectThemeColor = msoThemeColorAccent1 
.ForeColor.TintAndShade = 0 
.ForeColor. Brightness = -0.25 
.Transparency = 0 
-Visible = msoTrue 
.DashStyle = msoLineSysDot 

End With 

End Sub 


Figure A22. (Continued.) 
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